#' ---
#' title: "Replication: Can a Constitutional Monarch Influence Democratic Preferences? Japanese Emperor and the Regulation of Public Expression"
#' author: "Gento Kato and Susumu Annaka"
#' date: "March 16, 2022"
#' output:
#'  pdf_document: 
#'   latex_engine: xelatex 
#' toc: true
#' toc_depth: 4
#' ---

#+ echo=FALSE
# This document is encoded in UTF-8

#+ echo=FALSE
## Set Working Directory to the current directory 
## (If using RStudio, can be set automatically)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Clear Workspace
rm(list=ls())

## Used Packages
library(readr);library(knitr);library(ggplot2)
require(grid);library(gridExtra);library(haven)
require(sandwich); require(lmtest);require(texreg)

#+ message=FALSE
## Read Data
d <- readRDS("main_data_emperor_v2.rds")
dpanel <- readRDS("main_data_emperor_panel_v2.rds")
nrow(d)

dtmp <- d

#'
#' # Descriptive Statistics
#'
#'
#' ## DV: Regulating Expression in Public (Figure D.1)
#' 

ptab <- data.frame(pr = c(table(d$limitexp1[d$frame_left==1])/sum(table(d$limitexp1[d$frame_left==1])),
                          table(d$limitexp2[d$frame_left==1])/sum(table(d$limitexp2[d$frame_left==1])),
                          table(d$limitexp1[d$frame_right==1])/sum(table(d$limitexp1[d$frame_right==1])),
                          table(d$limitexp2[d$frame_right==1])/sum(table(d$limitexp2[d$frame_right==1]))),
                   v = rep(rev(c("Should regulate\n actively",
                                 "Should regulate\n if necessary",
                                 "Need careful\n  consideration",
                                 "Should not regulate\n as much as possible",
                                 "Should not regulate\n in any case")),2),
                   q = rep(c("After frame","After negative endorsement"), each=5),
                   f = rep(c("Hate Speech","Biased History"), each=10))
ptab$v <- factor(ptab$v, levels=unique(ptab$v))
ptab$q <- factor(ptab$q, levels=unique(ptab$q))
ptab$f <- factor(ptab$f, levels=unique(ptab$f))

p <- ggplot(ptab) + 
  geom_col(aes(x=v,y=pr,alpha=q), position = position_dodge(width=-1)) + 
  facet_grid(.~f) + 
  scale_alpha_manual("", values = c(0.4,0.8)) + 
  labs(y="Proportion",x="Response") +
  coord_flip() + 
  theme_bw() + 
  theme(legend.position = "bottom")

#+ fig.width=7, fig.height=5
p

#+ eval=FALSE
ggsave("../out/dvdist_v2.png", p, width=7, height=5)

#'
#' ## Ideologies (Figure E.2)
#'
#' 
#' ### Self-reported
#'
tab <- table(dtmp$ide_self)/sum(table(dtmp$ide_self))
tab <- data.frame(prop = as.numeric(tab),
                  names = c("Left\n(-3)\n","-2","-1",
                            "Mod.\n(0)\n","1","2","Right\n(3)\n"))
tab$names <- factor(tab$names, levels=tab$names)

p1 <- ggplot(tab, aes(x=names,y=prop)) + 
  geom_bar(stat="identity") + 
  ylab(NULL) + xlab(NULL) + 
  ggtitle("Self-reported\n(Frequency)") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p1

#' 
#' ### National security
#' 

p2_1 <- ggplot(dtmp[,"ide_iss_1"],
               aes(x=ide_iss_1,y=..count../sum(..count..))) +
  geom_histogram(bins=10, color="white") +
  ylab(NULL) + xlab(NULL) + 
  ggtitle("Nat. Security\n(Histogram)") + 
  scale_x_continuous(breaks=c(-3,-2,-1,0,1,2,3),
                     limits=c(-3,3),
                     labels=c("Left\n(-3)\n","-2","-1","0","1","2","Right\n(3)\n")) +
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p2_1

#' 
#' ### Equality
#' 

p2_2 <- ggplot(dtmp[,"ide_iss_2"],
               aes(x=ide_iss_2,y=..count../sum(..count..))) + 
  geom_histogram(bins=10,color="white") +
  ylab(NULL) + xlab(NULL) + 
  ggtitle("Equality\n(Histogram)") + 
  scale_x_continuous(breaks=c(-3,-2,-1,0,1,2,3),
                     limits=c(-3,3),
                     labels=c("Left\n(-3)\n","-2","-1","0","1","2","Right\n(3)\n")) +
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p2_2

#'
#' ### Party Support
#' 
tab <- table(dtmp$ide_psup)/sum(table(dtmp$ide_psup))
tab <- data.frame(prop = as.numeric(tab),
                  names = c("Left\nParty\n(-1)","No\nParty\n(0)",
                            "Right\nParty\n(1)"))
tab$names <- factor(tab$names, levels=tab$names)

p3 <- ggplot(tab, aes(x=names,y=prop)) + 
  geom_bar(stat="identity") + 
  ylab(NULL) + xlab(NULL) + 
  ggtitle("Party\n(Frequency)") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p3

#+ fig.width=15, fig.height=6
ggplot() + theme_void()
p <- arrangeGrob(p1,p2_1,p2_2,p3, nrow=1, left="Proportion")
grid.draw(p)

#+ eval=FALSE
ggsave("../out/idedist123_v2.png", p, width=15, height=6)

#'
#' # Analysis
#' 
#' ## Prepare Variable List
#' 

vn <- list("(Intercept)" = "(Intercept)",
           "frame_left" = "Hate speech frame",
           "frame_right" = "Biased history frame",
           "after" = "After endorsement",
           "Aword_expert" = "Neg. endorsement (scholars)",
           "Aword_emperor" = "Neg. endorsement (emperor)",
           "Aword_emperor_liberal" = "Neg. endorsement (liberal emperor)",
           "Bword_expert" = "Scholar (before treatment)",
           "Bword_emperor" = "Emperor (before endorsement)",
           "Bword_emperor_liberal" = "Liberal emperor (before endorsement)",
           "word_expert" = "Neg. endorsement (scholars)",
           "word_emperor" = "Neg. endorsement (emperor)",
           "word_emperor_liberal" = "Neg. endorsement (liberal emperor)",
           "frame_left:ide_self" = "Hate speech * ideology (self-reported)",
           "frame_left:ide_iss_1" = "Hate speech * ideology (national security)",
           "frame_left:ide_iss_2" = "Hate speech * ideology (equality)",
           "frame_left:left_psup" = "Hate speech * left party",
           "frame_left:right_psup" = "Hate speech * right party",
           "frame_right:ide_self" = "Biased history * ideology (self-reported)",
           "frame_right:ide_iss_1" = "Biased history * ideology (national security)",
           "frame_right:ide_iss_2" = "Biased history * ideology (equality)",
           "frame_right:left_psup" = "Biased history * left party",
           "frame_right:right_psup" = "Biased history * right party",
           "ide_self" = "Ideology (self-reported)",
           "ide_iss_1" = "Ideology (national security)",
           "ide_iss_2" = "Ideology (equality)",
           "left_psup" = "Left party support",
           "right_psup" = "Right party support",
           "fem" = "Gender (female)",
           "age" = "Age",
           "inccatMiddle (>=4m,<8m)" = "Income (middle)",
           "inccatHigh (>=8m)" = "Income (high)",
           "inccatMissing" = "Income (missing)",
           "educat>SHS & <University" = "Education (junior college/tech. school)",
           "educat>=University" = "Education (university)",
           "employed" = "Employed",
           "knall" = "Political knowledge (0-1)",
           "csup" = "Approve Abe Cabinet",
           "eqview" = "Japanese society is equal (0-1)",
           "hmview" = "Japanese society is homogeneous (0-1)",
           "0|1" = "Cut: No, any case|No, as much as possible",
           "1|2" = "Cut: No, as much as possible|Need to be careful",
           "2|3" = "Cut: Need to be careful|Yes, if necessary",
           "3|4" = "Cut: Yes, if necessary|Yes, actively",
           "-2|-1" = "Cut: -2 or under|-1",
           "-1|0" = "Cut: -1|0", 
           "0|2" = "Cut: 0|1 or above",
           "Weaker|Same" = "Cut: Weaker|Same",
           "Same|Stronger" = "Cut: Same|Stronger")

require(MASS)

#'
#' ## Check Conditional Framing Effect
#'
#' ### T-test (Table F.1)
#'

## Self-reported Ideology
table(d$ide_self)
mtf01_L <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_self>0)
mtf01_M <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_self==0)
mtf01_R <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_self<0)

### Issue Ideology 1 (National Security)
summary(d$ide_iss_1); sd(d$ide_iss_1)
mtf02A_L <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_iss_1> 1)
mtf02A_M <- t.test(limitexp1 ~ frame_right, data=d, subset=abs(ide_iss_1)<= 1)
mtf02A_R <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_iss_1< -1)

### Issue Ideology 2 (Equality)
summary(d$ide_iss_2); sd(d$ide_iss_2)
mtf02B_L <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_iss_2> 1)
mtf02B_M <- t.test(limitexp1 ~ frame_right, data=d, subset=abs(ide_iss_2)<= 1)
mtf02B_R <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_iss_2< -1)

### Party Support
table(d$ide_psup)
mtf03_L <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_psup==1)
mtf03_M <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_psup==0)
mtf03_R <- t.test(limitexp1 ~ frame_right, data=d, subset=ide_psup==-1)


mtf <- data.frame(Measurement = c("Self-reported","","", 
                           "Nat. Security","","", 
                           "Equality","","", 
                           "Party Support","",""),
           Ideology = c("Right (>0)","Moderate (0)","Left (<0)",
                        "Right (>1)","Moderate (-1:1)","Left (<-1)",
                        "Right (>1)","Moderate (-1:1)","Left (<-1)",
                        "Right (1)","Moderate (0)","Left (-1)"),
           "Hate Speech Mean" = 
             sprintf("%.3f",c(mtf01_L$estimate[1],mtf01_M$estimate[1],mtf01_R$estimate[1],
             mtf02A_L$estimate[1],mtf02A_M$estimate[1],mtf02A_R$estimate[1],
             mtf02B_L$estimate[1],mtf02B_M$estimate[1],mtf02B_R$estimate[1],
             mtf03_L$estimate[1],mtf03_M$estimate[1],mtf03_R$estimate[1])),
           "Biased History Mean" = 
             sprintf("%.3f",c(mtf01_L$estimate[2],mtf01_M$estimate[2],mtf01_R$estimate[2],
               mtf02A_L$estimate[2],mtf02A_M$estimate[2],mtf02A_R$estimate[2],
               mtf02B_L$estimate[2],mtf02B_M$estimate[2],mtf02B_R$estimate[2],
               mtf03_L$estimate[2],mtf03_M$estimate[2],mtf03_R$estimate[2])),
           p = sprintf("%.3f",c(mtf01_L$p.value,mtf01_M$p.value,mtf01_R$p.value,
                 mtf02A_L$p.value,mtf02A_M$p.value,mtf02A_R$p.value,
                 mtf02B_L$p.value,mtf02B_M$p.value,mtf02B_R$p.value,
                 mtf03_L$p.value,mtf03_M$p.value,mtf03_R$p.value)), 
           check.names = FALSE)

require(stargazer)
stargazer(mtf, summary=F, align=T, rownames = F, 
          title = "The ideology-moderated framing treatment effect on support for regulating expression in public places (t-test)",
          type = "text")

#+ eval=FALSE
stargazer(mtf, summary=F, align=T, rownames = F, 
          title = "The ideology-moderated framing treatment effect on support for regulating expression in public places (t-test)", 
          out="../out/resout_frame_ttest_v2.tex", 
          type = "latex")

#'
#' ### OLS: Estimation
#' 
#' #### Baseline (Table F.2)
#' 

## Self-reported Ideology
mbfL01 <- lm(limitexp1 ~ 
               frame_right * ide_self, 
             data=d)
summary(mbfL01)
coeftest(mbfL01, vcov.=vcovCL(mbfL01,cluster=na.omit(d[,c("start_id",all.vars(mbfL01$terms))])$start_id))
mbfL01c_se <- sqrt(diag(vcovCL(mbfL01,cluster=na.omit(d[,c("start_id",all.vars(mbfL01$terms))])$start_id)))
mbfL01c_p <- pt(-abs(summary(mbfL01)$coefficients[,1]/mbfL01c_se), df = mbfL01$df.residual)*2

## Issue Ideology
mbfL02 <- lm(limitexp1 ~ 
               frame_right * ide_iss_1 + 
               frame_right * ide_iss_2,  
             data=d)
summary(mbfL02)
coeftest(mbfL02, vcov.=vcovCL(mbfL02,cluster=na.omit(d[,c("start_id",all.vars(mbfL02$terms))])$start_id))
mbfL02c_se <- sqrt(diag(vcovCL(mbfL02,cluster=na.omit(d[,c("start_id",all.vars(mbfL02$terms))])$start_id)))
mbfL02c_p <- pt(-abs(summary(mbfL02)$coefficients[,1]/mbfL02c_se), df = mbfL02$df.residual)*2

## Party Support
mbfL03 <- lm(limitexp1 ~ 
               frame_right * left_psup + 
               frame_right * right_psup, 
             data=d)
summary(mbfL03)
coeftest(mbfL03, vcov.=vcovCL(mbfL03,cluster=na.omit(d[,c("start_id",all.vars(mbfL03$terms))])$start_id))
mbfL03c_se <- sqrt(diag(vcovCL(mbfL03,cluster=na.omit(d[,c("start_id",all.vars(mbfL03$terms))])$start_id)))
mbfL03c_p <- pt(-abs(summary(mbfL03)$coefficients[,1]/mbfL03c_se), df = mbfL03$df.residual)*2

screenreg(list(mbfL01,mbfL02,mbfL03), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbfL01c_se,mbfL02c_se,mbfL03c_se),
          override.pvalues = list(mbfL01c_p,mbfL02c_p,mbfL03c_p),
          custom.model.names = c("Self-reported","Issue","Party"),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mbfL01,mbfL02,mbfL03), stars = c(0.1,0.05,0.01,0.001), symbol = "\\dagger", 
       beside = T, digits = 3, single.row = T,
       override.se = list(mbfL01c_se,mbfL02c_se,mbfL03c_se),
       override.pvalues = list(mbfL01c_p,mbfL02c_p,mbfL03c_p),
       custom.model.names = c("Self-reported","Issue","Party"),
       custom.coef.map = vn, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "The ideology-moderated framing treatment effect on support for regulating expression in public places (OLS, baseline models)",
       file = "../out/resout_frame_lm_base_v2.tex", label = "table:resout_frame_lm_base")

#'
#' #### Extended (Table F.3)
#'

## Self-reported Ideology
mfL01 <- lm(limitexp1 ~ 
              frame_right * ide_self + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
            data=d)
summary(mfL01)
coeftest(mfL01, vcov.=vcovCL(mfL01,cluster=na.omit(d[,c("start_id",all.vars(mfL01$terms))])$start_id))
mfL01c_se <- sqrt(diag(vcovCL(mfL01,cluster=na.omit(d[,c("start_id",all.vars(mfL01$terms))])$start_id)))
mfL01c_p <- pt(-abs(summary(mfL01)$coefficients[,1]/mfL01c_se), df = mfL01$df.residual)*2

## Issue Ideology
mfL02 <- lm(limitexp1 ~ 
              frame_right * ide_iss_1 + 
              frame_right * ide_iss_2 + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
            data=d)
summary(mfL02)
coeftest(mfL02, vcov.=vcovCL(mfL02,cluster=na.omit(d[,c("start_id",all.vars(mfL02$terms))])$start_id))
mfL02c_se <- sqrt(diag(vcovCL(mfL02,cluster=na.omit(d[,c("start_id",all.vars(mfL02$terms))])$start_id)))
mfL02c_p <- pt(-abs(summary(mfL02)$coefficients[,1]/mfL02c_se), df = mfL02$df.residual)*2

## Party Support
mfL03 <- lm(limitexp1 ~ 
              frame_right * left_psup + 
              frame_right * right_psup + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
            data=d)
summary(mfL03)
coeftest(mfL03, vcov.=vcovCL(mfL03,cluster=na.omit(d[,c("start_id",all.vars(mfL03$terms))])$start_id))
mfL03c_se <- sqrt(diag(vcovCL(mfL03,cluster=na.omit(d[,c("start_id",all.vars(mfL03$terms))])$start_id)))
mfL03c_p <- pt(-abs(summary(mfL03)$coefficients[,1]/mfL03c_se), df = mfL03$df.residual)*2

screenreg(list(mfL01,mfL02,mfL03), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mfL01c_se,mfL02c_se,mfL03c_se),
          override.pvalues = list(mfL01c_p,mfL02c_p,mfL03c_p),
          custom.model.names = c("Self-reported","Issue","Party"),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mfL01,mfL02,mfL03), stars = c(0.1,0.05,0.01,0.001), symbol = "\\dagger", 
       beside = T, digits = 3, single.row = T,
       override.se = list(mfL01c_se,mfL02c_se,mfL03c_se),
       override.pvalues = list(mfL01c_p,mfL02c_p,mfL03c_p),
       custom.model.names = c("Self-reported","Issue","Party"),
       custom.coef.map = vn, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "The ideology-moderated framing treatment effect on support for regulating expression in public places (OLS, extended models)",
       file = "../out/resout_frame_lm_v2.tex", label = "table:resout_frame_lm")

#'
#' ### OLS: Conditional Effect Sizes
#' 
#' #### Baseline Models (Figure 1)
#' 

quantile(d$ide_self, probs = c(0.05,0.5,0.95))

## Self-reported Ideology ##

tmp1 <- lm(limitexp1 ~ 
             frame_right * I(ide_self+2), 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * I(ide_self-2),
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mbfL01$coefficients[2],mbfL01c_se[2],
          coefci(mbfL01, vcov.=vcovCL(mbfL01,cluster=na.omit(d[,c("start_id",all.vars(mbfL01$terms))])$start_id), level=0.95)[2,],
          coefci(mbfL01, vcov.=vcovCL(mbfL01,cluster=na.omit(d[,c("start_id",all.vars(mbfL01$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfL01dt <- tmpdt
mbfL01dt$type <- "Self-reported"
mbfL01dt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbfL01dt

## Issue Ideology (National Security) ##

quantile(d$ide_iss_1, probs = c(0.05,0.5,0.95))

tmp1 <- lm(limitexp1 ~ 
             frame_right * I(ide_iss_1+1.9) + 
             frame_right * ide_iss_2, 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * I(ide_iss_1-2) + 
             frame_right * ide_iss_2, 
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mbfL02$coefficients[2],mbfL02c_se[2],
          coefci(mbfL02, vcov.=vcovCL(mbfL02,cluster=na.omit(d[,c("start_id",all.vars(mbfL02$terms))])$start_id), level=0.95)[2,],
          coefci(mbfL02, vcov.=vcovCL(mbfL02,cluster=na.omit(d[,c("start_id",all.vars(mbfL02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfL02adt <- tmpdt
mbfL02adt$type <- "National\nSecurity"
mbfL02adt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbfL02adt

## Issue Ideology (Equality) ##

quantile(d$ide_iss_2, probs = c(0.05,0.5,0.95))

tmp1 <- lm(limitexp1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2+1.7), 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2-2.1), 
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mbfL02$coefficients[2],mbfL02c_se[2],
          coefci(mbfL02, vcov.=vcovCL(mbfL02,cluster=na.omit(d[,c("start_id",all.vars(mbfL02$terms))])$start_id), level=0.95)[2,],
          coefci(mbfL02, vcov.=vcovCL(mbfL02,cluster=na.omit(d[,c("start_id",all.vars(mbfL02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfL02bdt <- tmpdt
mbfL02bdt$type <- "Equality"
mbfL02bdt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbfL02bdt

## Party Support Ideology ##

tmp1 <- lm(limitexp1 ~ 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
             frame_right * right_psup, 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * left_psup + 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0), 
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mbfL03$coefficients[2],mbfL03c_se[2],
          coefci(mbfL03, vcov.=vcovCL(mbfL03,cluster=na.omit(d[,c("start_id",all.vars(mbfL03$terms))])$start_id), level=0.95)[2,],
          coefci(mbfL03, vcov.=vcovCL(mbfL03,cluster=na.omit(d[,c("start_id",all.vars(mbfL03$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfL03dt <- tmpdt
mbfL03dt$type <- "Party"
mbfL03dt$value <- c("Neither","Left\nParty","Right\nParty")
mbfL03dt

## Combine Data ##

mbfLdt <- rbind(mbfL01dt,mbfL02adt,mbfL02bdt,mbfL03dt)
mbfLdt$type <- factor(mbfLdt$type, unique(mbfLdt$type))
unique(mbfLdt$value)
mbfLdt$value <- factor(mbfLdt$value, c("Left\n(5%)", "Left\nParty",
                                       "Neut.\n(50%)", "Neither",
                                       "Right\n(95%)", "Right\nParty"))

write.csv(mbfLdt, row.names = F, file = "../out/effect_frame_lm_base_v2.csv")

## Plot of Conditional Framing Effect ##

require(ggplot2)

p <- ggplot(mbfLdt, aes(x=value,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~type, scales = "free_x") + 
  labs(x="Ideological position", y="Biased history vs. hate speech frame\nOLS coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_frame_lm_base_v2.pdf", width=6, height=4)
ggsave("../out/effect_frame_lm_base_v2.png", width=6, height=4)

#' 
#' #### Extended Models (Figure F.1)
#' 

quantile(d$ide_self, probs = c(0.05,0.5,0.95))

## Self-reported Ideology ##

tmp1 <- lm(limitexp1 ~ 
             frame_right * I(ide_self+2) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * I(ide_self-2) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mfL01$coefficients[2],mfL01c_se[2],
          coefci(mfL01, vcov.=vcovCL(mfL01,cluster=na.omit(d[,c("start_id",all.vars(mfL01$terms))])$start_id), level=0.95)[2,],
          coefci(mfL01, vcov.=vcovCL(mfL01,cluster=na.omit(d[,c("start_id",all.vars(mfL01$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfL01dt <- tmpdt
mfL01dt$type <- "Self-reported"
mfL01dt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mfL01dt

## Issue Ideology (National Security) ##

quantile(d$ide_iss_1, probs = c(0.05,0.5,0.95))

tmp1 <- lm(limitexp1 ~ 
             frame_right * I(ide_iss_1+1.9) + 
             frame_right * ide_iss_2 + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * I(ide_iss_1-2) + 
             frame_right * ide_iss_2 + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mfL02$coefficients[2],mfL02c_se[2],
          coefci(mfL02, vcov.=vcovCL(mfL02,cluster=na.omit(d[,c("start_id",all.vars(mfL02$terms))])$start_id), level=0.95)[2,],
          coefci(mfL02, vcov.=vcovCL(mfL02,cluster=na.omit(d[,c("start_id",all.vars(mfL02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfL02adt <- tmpdt
mfL02adt$type <- "National\nSecurity"
mfL02adt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mfL02adt

## Issue Ideology (Equality) ##

quantile(d$ide_iss_2, probs = c(0.05,0.5,0.95))

tmp1 <- lm(limitexp1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2+1.7) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2-2.1) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mfL02$coefficients[2],mfL02c_se[2],
          coefci(mfL02, vcov.=vcovCL(mfL02,cluster=na.omit(d[,c("start_id",all.vars(mfL02$terms))])$start_id), level=0.95)[2,],
          coefci(mfL02, vcov.=vcovCL(mfL02,cluster=na.omit(d[,c("start_id",all.vars(mfL02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfL02bdt <- tmpdt
mfL02bdt$type <- "Equality"
mfL02bdt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mfL02bdt

## Party Support Ideology ##

tmp1 <- lm(limitexp1 ~ 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
             frame_right * right_psup + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp1 ~ 
             frame_right * left_psup + 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt <- 
  rbind(c(mfL03$coefficients[2],mfL03c_se[2],
          coefci(mfL03, vcov.=vcovCL(mfL03,cluster=na.omit(d[,c("start_id",all.vars(mfL03$terms))])$start_id), level=0.95)[2,],
          coefci(mfL03, vcov.=vcovCL(mfL03,cluster=na.omit(d[,c("start_id",all.vars(mfL03$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfL03dt <- tmpdt
mfL03dt$type <- "Party"
mfL03dt$value <- c("Neither","Left\nParty","Right\nParty")
mfL03dt

## Combine Data ##

mfLdt <- rbind(mfL01dt,mfL02adt,mfL02bdt,mfL03dt)
mfLdt$type <- factor(mfLdt$type, unique(mfLdt$type))
unique(mfLdt$value)
mfLdt$value <- factor(mfLdt$value, c("Left\n(5%)", "Left\nParty",
                                       "Neut.\n(50%)", "Neither",
                                       "Right\n(95%)", "Right\nParty"))

write.csv(mfLdt, row.names = F, file = "../out/effect_frame_lm_v2.csv")

## Plot of Conditional Framing Effect ##

require(ggplot2)

p <- ggplot(mfLdt, aes(x=value,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~type, scales = "free_x") + 
  labs(x="Ideological position", y="Biased history vs. hate speech frame\nOLS coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_frame_lm_v2.pdf", width=6, height=4)
ggsave("../out/effect_frame_lm_v2.png", width=6, height=4)

#'
#' ### Logit: Estimation
#' 
#' #### Baseline (Table F.4)
#' 

## Self-reported Ideology
mbfG01 <- glm(limitexp1>1 ~ 
               frame_right * ide_self, 
             data=d, family=binomial("logit"))
summary(mbfG01)
coeftest(mbfG01, vcov.=vcovCL(mbfG01,cluster=na.omit(d[,c("start_id",all.vars(mbfG01$terms))])$start_id))
mbfG01c_se <- sqrt(diag(vcovCL(mbfG01,cluster=na.omit(d[,c("start_id",all.vars(mbfG01$terms))])$start_id)))
mbfG01c_p <- pnorm(-abs(summary(mbfG01)$coefficients[,1]/mbfG01c_se))*2

## Issue Ideology
mbfG02 <- glm(limitexp1>1 ~ 
               frame_right * ide_iss_1 + 
               frame_right * ide_iss_2,  
             data=d, family=binomial("logit"))
summary(mbfG02)
coeftest(mbfG02, vcov.=vcovCL(mbfG02,cluster=na.omit(d[,c("start_id",all.vars(mbfG02$terms))])$start_id))
mbfG02c_se <- sqrt(diag(vcovCL(mbfG02,cluster=na.omit(d[,c("start_id",all.vars(mbfG02$terms))])$start_id)))
mbfG02c_p <- pnorm(-abs(summary(mbfG02)$coefficients[,1]/mbfG02c_se))*2

## Party Support
mbfG03 <- glm(limitexp1>1 ~ 
               frame_right * left_psup + 
               frame_right * right_psup, 
             data=d, family=binomial("logit"))
summary(mbfG03)
coeftest(mbfG03, vcov.=vcovCL(mbfG03,cluster=na.omit(d[,c("start_id",all.vars(mbfG03$terms))])$start_id))
mbfG03c_se <- sqrt(diag(vcovCL(mbfG03,cluster=na.omit(d[,c("start_id",all.vars(mbfG03$terms))])$start_id)))
mbfG03c_p <- pnorm(-abs(summary(mbfG03)$coefficients[,1]/mbfG03c_se))*2

screenreg(list(mbfG01,mbfG02,mbfG03), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbfG01c_se,mbfG02c_se,mbfG03c_se),
          override.pvalues = list(mbfG01c_p,mbfG02c_p,mbfG03c_p),
          custom.model.names = c("Self-reported","Issue","Party"),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mbfG01,mbfG02,mbfG03), stars = c(0.1,0.05,0.01,0.001), symbol = "\\dagger", 
       beside = T, digits = 3, single.row = T,
       override.se = list(mbfG01c_se,mbfG02c_se,mbfG03c_se),
       override.pvalues = list(mbfG01c_p,mbfG02c_p,mbfG03c_p),
       custom.model.names = c("Self-reported","Issue","Party"),
       custom.coef.map = vn, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "The ideology-moderated framing treatment effect on support for regulating expression in public places (OLS, baseline models)",
       file = "../out/resout_frame_logit_base_v2.tex", label = "table:resout_frame_logit_base")

#'
#' #### Extended (Table F.5)
#'

## Self-reported Ideology
mfG01 <- glm(limitexp1>1 ~ 
              frame_right * ide_self + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
            data=d, family=binomial("logit"))
summary(mfG01)
coeftest(mfG01, vcov.=vcovCL(mfG01,cluster=na.omit(d[,c("start_id",all.vars(mfG01$terms))])$start_id))
mfG01c_se <- sqrt(diag(vcovCL(mfG01,cluster=na.omit(d[,c("start_id",all.vars(mfG01$terms))])$start_id)))
mfG01c_p <- pnorm(-abs(summary(mfG01)$coefficients[,1]/mfG01c_se))*2

## Issue Ideology
mfG02 <- glm(limitexp1>1 ~ 
              frame_right * ide_iss_1 + 
              frame_right * ide_iss_2 + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
            data=d, family=binomial("logit"))
summary(mfG02)
coeftest(mfG02, vcov.=vcovCL(mfG02,cluster=na.omit(d[,c("start_id",all.vars(mfG02$terms))])$start_id))
mfG02c_se <- sqrt(diag(vcovCL(mfG02,cluster=na.omit(d[,c("start_id",all.vars(mfG02$terms))])$start_id)))
mfG02c_p <- pnorm(-abs(summary(mfG02)$coefficients[,1]/mfG02c_se))*2

## Party Support
mfG03 <- glm(limitexp1>1 ~ 
              frame_right * left_psup + 
              frame_right * right_psup + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
            data=d, family=binomial("logit"))
summary(mfG03)
coeftest(mfG03, vcov.=vcovCL(mfG03,cluster=na.omit(d[,c("start_id",all.vars(mfG03$terms))])$start_id))
mfG03c_se <- sqrt(diag(vcovCL(mfG03,cluster=na.omit(d[,c("start_id",all.vars(mfG03$terms))])$start_id)))
mfG03c_p <- pnorm(-abs(summary(mfG03)$coefficients[,1]/mfG03c_se))*2

screenreg(list(mfG01,mfG02,mfG03), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mfG01c_se,mfG02c_se,mfG03c_se),
          override.pvalues = list(mfG01c_p,mfG02c_p,mfG03c_p),
          custom.model.names = c("Self-reported","Issue","Party"),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mfG01,mfG02,mfG03), stars = c(0.1,0.05,0.01,0.001), symbol = "\\dagger", 
       beside = T, digits = 3, single.row = T,
       override.se = list(mfG01c_se,mfG02c_se,mfG03c_se),
       override.pvalues = list(mfG01c_p,mfG02c_p,mfG03c_p),
       custom.model.names = c("Self-reported","Issue","Party"),
       custom.coef.map = vn, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "The ideology-moderated framing treatment effect on support for regulating expression in public places (OLS, extended models)",
       file = "../out/resout_frame_logit_v2.tex", label = "table:resout_frame_logit")

#'
#' ### Logit: Conditional Effect Sizes
#' 
#' #### Baseline Models (Figure F.2)
#' 

quantile(d$ide_self, probs = c(0.05,0.5,0.95))

## Self-reported Ideology ##

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_self+2), 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_self-2),
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbfG01$coefficients[2],mbfG01c_se[2],
          coefci(mbfG01, vcov.=vcovCL(mbfG01,cluster=na.omit(d[,c("start_id",all.vars(mbfG01$terms))])$start_id), level=0.95)[2,],
          coefci(mbfG01, vcov.=vcovCL(mbfG01,cluster=na.omit(d[,c("start_id",all.vars(mbfG01$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfG01dt <- tmpdt
mbfG01dt$type <- "Self-reported"
mbfG01dt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbfG01dt

## Issue Ideology (National Security) ##

quantile(d$ide_iss_1, probs = c(0.05,0.5,0.95))

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_iss_1+1.9) + 
             frame_right * ide_iss_2, 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_iss_1-2) + 
             frame_right * ide_iss_2, 
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbfG02$coefficients[2],mbfG02c_se[2],
          coefci(mbfG02, vcov.=vcovCL(mbfG02,cluster=na.omit(d[,c("start_id",all.vars(mbfG02$terms))])$start_id), level=0.95)[2,],
          coefci(mbfG02, vcov.=vcovCL(mbfG02,cluster=na.omit(d[,c("start_id",all.vars(mbfG02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfG02adt <- tmpdt
mbfG02adt$type <- "National\nSecurity"
mbfG02adt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbfG02adt

## Issue Ideology (Equality) ##

quantile(d$ide_iss_2, probs = c(0.05,0.5,0.95))

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2+1.7), 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2-2.1), 
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbfG02$coefficients[2],mbfG02c_se[2],
          coefci(mbfG02, vcov.=vcovCL(mbfG02,cluster=na.omit(d[,c("start_id",all.vars(mbfG02$terms))])$start_id), level=0.95)[2,],
          coefci(mbfG02, vcov.=vcovCL(mbfG02,cluster=na.omit(d[,c("start_id",all.vars(mbfG02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfG02bdt <- tmpdt
mbfG02bdt$type <- "Equality"
mbfG02bdt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbfG02bdt

## Party Support Ideology ##

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
             frame_right * right_psup, 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * left_psup + 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0), 
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbfG03$coefficients[2],mbfG03c_se[2],
          coefci(mbfG03, vcov.=vcovCL(mbfG03,cluster=na.omit(d[,c("start_id",all.vars(mbfG03$terms))])$start_id), level=0.95)[2,],
          coefci(mbfG03, vcov.=vcovCL(mbfG03,cluster=na.omit(d[,c("start_id",all.vars(mbfG03$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbfG03dt <- tmpdt
mbfG03dt$type <- "Party"
mbfG03dt$value <- c("Neither","Left\nParty","Right\nParty")
mbfG03dt

## Combine Data ##

mbfGdt <- rbind(mbfG01dt,mbfG02adt,mbfG02bdt,mbfG03dt)
mbfGdt$type <- factor(mbfGdt$type, unique(mbfGdt$type))
unique(mbfGdt$value)
mbfGdt$value <- factor(mbfGdt$value, c("Left\n(5%)", "Left\nParty",
                                       "Neut.\n(50%)", "Neither",
                                       "Right\n(95%)", "Right\nParty"))

write.csv(mbfGdt, row.names = F, file = "../out/effect_frame_logit_base_v2.csv")

## Plot of Conditional Framing Effect ##

require(ggplot2)

p <- ggplot(mbfGdt, aes(x=value,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~type, scales = "free_x") + 
  labs(x="Ideological position", y="Biased history vs. hate speech frame\nOLS coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_frame_logit_base_v2.pdf", width=6, height=4)
ggsave("../out/effect_frame_logit_base_v2.png", width=6, height=4)

#' 
#' #### Extended Models (Figure F.3)
#' 

quantile(d$ide_self, probs = c(0.05,0.5,0.95))

## Self-reported Ideology ##

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_self+2) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_self-2) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mfG01$coefficients[2],mfG01c_se[2],
          coefci(mfG01, vcov.=vcovCL(mfG01,cluster=na.omit(d[,c("start_id",all.vars(mfG01$terms))])$start_id), level=0.95)[2,],
          coefci(mfG01, vcov.=vcovCL(mfG01,cluster=na.omit(d[,c("start_id",all.vars(mfG01$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfG01dt <- tmpdt
mfG01dt$type <- "Self-reported"
mfG01dt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mfG01dt

## Issue Ideology (National Security) ##

quantile(d$ide_iss_1, probs = c(0.05,0.5,0.95))

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_iss_1+1.9) + 
             frame_right * ide_iss_2 + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * I(ide_iss_1-2) + 
             frame_right * ide_iss_2 + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mfG02$coefficients[2],mfG02c_se[2],
          coefci(mfG02, vcov.=vcovCL(mfG02,cluster=na.omit(d[,c("start_id",all.vars(mfG02$terms))])$start_id), level=0.95)[2,],
          coefci(mfG02, vcov.=vcovCL(mfG02,cluster=na.omit(d[,c("start_id",all.vars(mfG02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfG02adt <- tmpdt
mfG02adt$type <- "National\nSecurity"
mfG02adt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mfG02adt

## Issue Ideology (Equality) ##

quantile(d$ide_iss_2, probs = c(0.05,0.5,0.95))

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2+1.7) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * ide_iss_1 + 
             frame_right * I(ide_iss_2-2.1) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mfG02$coefficients[2],mfG02c_se[2],
          coefci(mfG02, vcov.=vcovCL(mfG02,cluster=na.omit(d[,c("start_id",all.vars(mfG02$terms))])$start_id), level=0.95)[2,],
          coefci(mfG02, vcov.=vcovCL(mfG02,cluster=na.omit(d[,c("start_id",all.vars(mfG02$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfG02bdt <- tmpdt
mfG02bdt$type <- "Equality"
mfG02bdt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mfG02bdt

## Party Support Ideology ##

tmp1 <- glm(limitexp1>1 ~ 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
             frame_right * right_psup + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp1>1 ~ 
             frame_right * left_psup + 
             frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data=d, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mfG03$coefficients[2],mfG03c_se[2],
          coefci(mfG03, vcov.=vcovCL(mfG03,cluster=na.omit(d[,c("start_id",all.vars(mfG03$terms))])$start_id), level=0.95)[2,],
          coefci(mfG03, vcov.=vcovCL(mfG03,cluster=na.omit(d[,c("start_id",all.vars(mfG03$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mfG03dt <- tmpdt
mfG03dt$type <- "Party"
mfG03dt$value <- c("Neither","Left\nParty","Right\nParty")
mfG03dt

## Combine Data ##

mfGdt <- rbind(mfG01dt,mfG02adt,mfG02bdt,mfG03dt)
mfGdt$type <- factor(mfGdt$type, unique(mfGdt$type))
unique(mfGdt$value)
mfGdt$value <- factor(mfGdt$value, c("Left\n(5%)", "Left\nParty",
                                     "Neut.\n(50%)", "Neither",
                                     "Right\n(95%)", "Right\nParty"))

write.csv(mfGdt, row.names = F, file = "../out/effect_frame_logit_v2.csv")

## Plot of Conditional Framing Effect ##

require(ggplot2)

p <- ggplot(mfGdt, aes(x=value,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~type, scales = "free_x") + 
  labs(x="Ideological position", y="Biased history vs. hate speech frame\nOLS coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_frame_logit_v2.pdf", width=6, height=4)
ggsave("../out/effect_frame_logit_v2.png", width=6, height=4)

#'
#' ### Ordinal Logit: Estimation
#' 
#' #### Baseline (Table F.6)
#' 

## Self-reported Ideology
mbf01 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ide_self, data=d, Hess = TRUE)
summary(mbf01)
coeftest(mbf01, vcov.=vcovCL(mbf01,cluster=na.omit(d[,c("start_id",all.vars(mbf01$terms))])$start_id))
mbf01c_se <- sqrt(diag(vcovCL(mbf01,cluster=na.omit(d[,c("start_id",all.vars(mbf01$terms))])$start_id)))
mbf01c_p <- pnorm(-abs(summary(mbf01)$coefficients[,1]/mbf01c_se))*2

## Issue Ideology
mbf02 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ide_iss_1 + 
               frame_right * ide_iss_2, data=d, Hess = TRUE)
summary(mbf02)
coeftest(mbf02, vcov.=vcovCL(mbf02,cluster=na.omit(d[,c("start_id",all.vars(mbf02$terms))])$start_id))
mbf02c_se <- sqrt(diag(vcovCL(mbf02,cluster=na.omit(d[,c("start_id",all.vars(mbf02$terms))])$start_id)))
mbf02c_p <- pnorm(-abs(summary(mbf02)$coefficients[,1]/mbf02c_se))*2

## Party Support
mbf03 <- polr(as.ordered(limitexp1) ~ 
               frame_right * left_psup + 
               frame_right * right_psup, data=d, Hess = TRUE)
summary(mbf03)
coeftest(mbf03, vcov.=vcovCL(mbf03,cluster=na.omit(d[,c("start_id",all.vars(mbf03$terms))])$start_id))
mbf03c_se <- sqrt(diag(vcovCL(mbf03,cluster=na.omit(d[,c("start_id",all.vars(mbf03$terms))])$start_id)))
mbf03c_p <- pnorm(-abs(summary(mbf03)$coefficients[,1]/mbf03c_se))*2

screenreg(list(mbf01,mbf02,mbf03), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbf01c_se,mbf02c_se,mbf03c_se),
          override.pvalues = list(mbf01c_p,mbf02c_p,mbf03c_p),
          custom.model.names = c("Self-reported","Issue","Party"),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mbf01,mbf02,mbf03), stars = c(0.1,0.05,0.01,0.001), symbol = "\\dagger", 
       beside = T, digits = 3, single.row = T,
       override.se = list(mbf01c_se,mbf02c_se,mbf03c_se),
       override.pvalues = list(mbf01c_p,mbf02c_p,mbf03c_p),
       custom.model.names = c("Self-reported","Issue","Party"),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "The ideology-moderated framing treatment effect on support for regulating expression in public places (ordinal logit, baseline models)",
       file = "../out/resout_frame_ol_base_v2.tex", label = "table:resout_frame_ol_base")

#'
#' #### Extended (Table F.7)
#'

## Self-reported Ideology
mf01 <- polr(as.ordered(limitexp1) ~ 
            frame_right * ide_self + 
            fem + age + inccat + educat + employed + knall + csup + eqview + hmview, data=d, Hess = TRUE)
summary(mf01)
coeftest(mf01, vcov.=vcovCL(mf01,cluster=na.omit(d[,c("start_id",all.vars(mf01$terms))])$start_id))
mf01c_se <- sqrt(diag(vcovCL(mf01,cluster=na.omit(d[,c("start_id",all.vars(mf01$terms))])$start_id)))
mf01c_p <- pnorm(-abs(summary(mf01)$coefficients[,1]/mf01c_se))*2

## Issue Ideology
mf02 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ide_iss_1 + 
               frame_right * ide_iss_2 + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, data=d, Hess = TRUE)
summary(mf02)
coeftest(mf02, vcov.=vcovCL(mf02,cluster=na.omit(d[,c("start_id",all.vars(mf02$terms))])$start_id))
mf02c_se <- sqrt(diag(vcovCL(mf02,cluster=na.omit(d[,c("start_id",all.vars(mf02$terms))])$start_id)))
mf02c_p <- pnorm(-abs(summary(mf02)$coefficients[,1]/mf02c_se))*2

## Party Support
mf03 <- polr(as.ordered(limitexp1) ~ 
               frame_right * left_psup + 
               frame_right * right_psup + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, data=d, Hess = TRUE)
summary(mf03)
coeftest(mf03, vcov.=vcovCL(mf03,cluster=na.omit(d[,c("start_id",all.vars(mf03$terms))])$start_id))
mf03c_se <- sqrt(diag(vcovCL(mf03,cluster=na.omit(d[,c("start_id",all.vars(mf03$terms))])$start_id)))
mf03c_p <- pnorm(-abs(summary(mf03)$coefficients[,1]/mf03c_se))*2

screenreg(list(mf01,mf02,mf03), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mf01c_se,mf02c_se,mf03c_se),
          override.pvalues = list(mf01c_p,mf02c_p,mf03c_p),
          custom.model.names = c("Self-reported","Issue","Party"),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mf01,mf02,mf03), stars = c(0.1,0.05,0.01,0.001), symbol = "\\dagger", 
       beside = T, digits = 3, single.row = T,
       override.se = list(mf01c_se,mf02c_se,mf03c_se),
       override.pvalues = list(mf01c_p,mf02c_p,mf03c_p),
       custom.model.names = c("Self-reported","Issue","Party"),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "The ideology-moderated framing treatment effect on support for regulating expression in public places (ordinal logit, extended models)",
       file = "../out/resout_frame_ol_v2.tex", label = "table:resout_frame_ol")

#'
#' ### Ordinal Logit: Conditional Effect Sizes
#'
#' #### Baseline Models (Figure F.4)
#'

quantile(d$ide_self, probs=c(0.05,0.5,0.95))

## Self-reported Ideology ##

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_self+2),
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_self-2), data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbf01$coefficients[1],mbf01c_se[1],
          coefci(mbf01, vcov.=vcovCL(mbf01,cluster=na.omit(d[,c("start_id",all.vars(mbf01$terms))])$start_id), level=0.95)[1,],
          coefci(mbf01, vcov.=vcovCL(mbf01,cluster=na.omit(d[,c("start_id",all.vars(mbf01$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbf01dt <- tmpdt
mbf01dt$type <- "Self-reported"
mbf01dt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbf01dt

## Issue Ideology (National Security) ##

quantile(d$ide_iss_1, probs = c(0.05, 0.5, 0.95))

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_iss_1+1.9) + 
               frame_right * ide_iss_2, 
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_iss_1-2) + 
               frame_right * ide_iss_2, 
             data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbf02$coefficients[1],mbf02c_se[1],
          coefci(mbf02, vcov.=vcovCL(mbf02,cluster=na.omit(d[,c("start_id",all.vars(mbf02$terms))])$start_id), level=0.95)[1,],
          coefci(mbf02, vcov.=vcovCL(mbf02,cluster=na.omit(d[,c("start_id",all.vars(mbf02$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbf02adt <- tmpdt
mbf02adt$type <- "National\nSecurity"
mbf02adt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbf02adt

## Issue Ideology (Equality) ##

quantile(d$ide_iss_2, probs = c(0.05,0.5,0.95))

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ide_iss_1 + 
               frame_right * I(ide_iss_2+1.7),
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ide_iss_1 + 
               frame_right * I(ide_iss_2-2.1),
             data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbf02$coefficients[1],mbf02c_se[1],
          coefci(mbf02, vcov.=vcovCL(mbf02,cluster=na.omit(d[,c("start_id",all.vars(mbf02$terms))])$start_id), level=0.95)[1,],
          coefci(mbf02, vcov.=vcovCL(mbf02,cluster=na.omit(d[,c("start_id",all.vars(mbf02$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbf02bdt <- tmpdt
mbf02bdt$type <- "Equality"
mbf02bdt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mbf02bdt

## Party Support Ideology ##

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
               frame_right * right_psup, 
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * left_psup + 
               frame_right * ifelse(left_psup==0&right_psup==0,1,0), 
             data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mbf03$coefficients[1],mbf03c_se[1],
          coefci(mbf03, vcov.=vcovCL(mbf03,cluster=na.omit(d[,c("start_id",all.vars(mbf03$terms))])$start_id), level=0.95)[1,],
          coefci(mbf03, vcov.=vcovCL(mbf03,cluster=na.omit(d[,c("start_id",all.vars(mbf03$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mbf03dt <- tmpdt
mbf03dt$type <- "Party"
mbf03dt$value <- c("Neither","Left\nParty","Right\nParty")
mbf03dt

## Combine Data ##

mbfdt <- rbind(mbf01dt,mbf02adt,mbf02bdt,mbf03dt)
mbfdt$type <- factor(mbfdt$type, unique(mbfdt$type))
unique(mbfdt$value)
mbfdt$value <- factor(mbfdt$value, c("Left\n(5%)", "Left\nParty",
                                   "Neut.\n(50%)", "Neither",
                                   "Right\n(95%)", "Right\nParty"))

write.csv(mbfdt, row.names = F, file = "../out/effect_frame_ol_base_v2.csv")

## Plot of Conditional Framing Effect ##

require(ggplot2)

p <- ggplot(mbfdt, aes(x=value,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~type, scales = "free_x") + 
  labs(x="Ideological position", y="Biased history vs. hate speech frame\nordinal logit coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_frame_ol_base_v2.pdf", width=6, height=4)
ggsave("../out/effect_frame_ol_base_v2.png", width=6, height=4)

#'
#' #### Extended Models (Figure F.5)
#'

quantile(d$ide_self, probs=c(0.05,0.5,0.95))

## Self-reported Ideology ##

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_self+2) + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_self-2) + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mf01$coefficients[1],mf01c_se[1],
          coefci(mf01, vcov.=vcovCL(mf01,cluster=na.omit(d[,c("start_id",all.vars(mf01$terms))])$start_id), level=0.95)[1,],
          coefci(mf01, vcov.=vcovCL(mf01,cluster=na.omit(d[,c("start_id",all.vars(mf01$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mf01dt <- tmpdt
mf01dt$type <- "Self-reported"
mf01dt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mf01dt

## Issue Ideology (National Security) ##

quantile(d$ide_iss_1, probs = c(0.05, 0.5, 0.95))

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_iss_1+1.9) + 
               frame_right * ide_iss_2 + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * I(ide_iss_1-2) + 
               frame_right * ide_iss_2 + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mf02$coefficients[1],mf02c_se[1],
          coefci(mf02, vcov.=vcovCL(mf02,cluster=na.omit(d[,c("start_id",all.vars(mf02$terms))])$start_id), level=0.95)[1,],
          coefci(mf02, vcov.=vcovCL(mf02,cluster=na.omit(d[,c("start_id",all.vars(mf02$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mf02adt <- tmpdt
mf02adt$type <- "National\nSecurity"
mf02adt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mf02adt

## Issue Ideology (Equality) ##

quantile(d$ide_iss_2, probs = c(0.05,0.5,0.95))

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ide_iss_1 + 
               frame_right * I(ide_iss_2+1.7) + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ide_iss_1 + 
               frame_right * I(ide_iss_2-2.1) + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview,
             data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mf02$coefficients[1],mf02c_se[1],
          coefci(mf02, vcov.=vcovCL(mf02,cluster=na.omit(d[,c("start_id",all.vars(mf02$terms))])$start_id), level=0.95)[1,],
          coefci(mf02, vcov.=vcovCL(mf02,cluster=na.omit(d[,c("start_id",all.vars(mf02$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mf02bdt <- tmpdt
mf02bdt$type <- "Equality"
mf02bdt$value <- c("Neut.\n(50%)","Left\n(5%)","Right\n(95%)")
mf02bdt

## Party Support Ideology ##

tmp1 <- polr(as.ordered(limitexp1) ~ 
               frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
               frame_right * right_psup + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data=d, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp1) ~ 
               frame_right * left_psup + 
               frame_right * ifelse(left_psup==0&right_psup==0,1,0) + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data=d, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt <- 
  rbind(c(mf03$coefficients[1],mf03c_se[1],
          coefci(mf03, vcov.=vcovCL(mf03,cluster=na.omit(d[,c("start_id",all.vars(mf03$terms))])$start_id), level=0.95)[1,],
          coefci(mf03, vcov.=vcovCL(mf03,cluster=na.omit(d[,c("start_id",all.vars(mf03$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(d[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(d[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt <- as.data.frame(tmpdt)
colnames(tmpdt) <- c("cf","se","lci","uci","lci90","uci90")

mf03dt <- tmpdt
mf03dt$type <- "Party"
mf03dt$value <- c("Neither","Left\nParty","Right\nParty")
mf03dt

## Combine Data ##

mfdt <- rbind(mf01dt,mf02adt,mf02bdt,mf03dt)
mfdt$type <- factor(mfdt$type, unique(mfdt$type))
unique(mfdt$value)
mfdt$value <- factor(mfdt$value, c("Left\n(5%)", "Left\nParty",
                                     "Neut.\n(50%)", "Neither",
                                     "Right\n(95%)", "Right\nParty"))

write.csv(mfdt, row.names = F, file = "../out/effect_frame_ol_v2.csv")

## Plot of Conditional Framing Effect ##

require(ggplot2)

p <- ggplot(mfdt, aes(x=value,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~type, scales = "free_x") + 
  labs(x="Ideological position", y="Biased history vs. hate speech frame\nordinal logit coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_frame_ol_v2.pdf", width=6, height=4)
ggsave("../out/effect_frame_ol_v2.png", width=6, height=4)

#'
#' ## Endorsement Treatment Effects 
#' 
#' ### Difference Models 
#'

## Subset data
dL <- subset(d, frame_left==1)
dR <- subset(d, frame_right==1)

#' 
#' #### T-test (Table G.1)
#'

## One Sample T-test
(mt0_L <- t.test(limitexpdif~1, data=dL, subset=word_expert==1))
(mt1_L <- t.test(limitexpdif~1, data=dL, subset=word_emperor==1))
(mt2_L <- t.test(limitexpdif~1, data=dL, subset=word_emperor_liberal==1))
(mt0_R <- t.test(limitexpdif~1, data=dR, subset=word_expert==1))
(mt1_R <- t.test(limitexpdif~1, data=dR, subset=word_emperor==1))
(mt2_R <- t.test(limitexpdif~1, data=dR, subset=word_emperor_liberal==1))

mt <- data.frame(Issue = c("Hate Speech","","","Biased History","",""),
                  Endorser = c("Experts","Emperor","Liberal Emperor"),
                  "Mean Difference" = 
                    sprintf("%.3f",c(mt0_L$estimate,mt1_L$estimate,mt2_L$estimate,
                                     mt0_R$estimate,mt1_R$estimate,mt2_R$estimate)),
                  p = sprintf("%.3f",c(mt0_L$p.value,mt1_L$p.value,mt2_L$p.value,
                                       mt0_R$p.value,mt1_R$p.value,mt2_R$p.value)), 
                  check.names = FALSE)

require(stargazer)
stargazer(mt, summary=F, align=T, rownames = F, 
          title = "Endorsement treatment effects on the difference in the support for regulating expression in public places (t-test)",
          type = "text")

#+ eval=FALSE
stargazer(mt, summary=F, align=T, rownames = F, 
          title = "Endorsement treatment effects on the difference in the support for regulating expression in public places (t-test)", 
          out="../out/endorsement_diff_ttest_v2.tex", 
          type = "latex")

#'
#' #### Linear Model (Baseline)
#'

## Baseline (Linear Model)
mbL01 <- lm(limitexpdif ~ word_emperor + word_emperor_liberal, data = dL)
coeftest(mbL01, vcov.=vcovCL(mbL01,cluster=na.omit(dL[,c("start_id",all.vars(mbL01$terms))])$start_id))
mbL01c_se <- sqrt(diag(vcovCL(mbL01,cluster=na.omit(dL[,c("start_id",all.vars(mbL01$terms))])$start_id)))
# coeftest(mbL01, vcov.=vcovHC(mbL01, type="HC2"))
# mbL01c_se <- sqrt(diag(vcovHC(mbL01, type="HC2")))
mbL01c_p <- pt(-abs(summary(mbL01)$coefficients[,1]/mbL01c_se), df = mbL01$df.residual)*2
mbL02 <- lm(limitexpdif ~ word_emperor + word_emperor_liberal, data = dR)
coeftest(mbL02, vcov.=vcovCL(mbL02,cluster=na.omit(dR[,c("start_id",all.vars(mbL02$terms))])$start_id))
mbL02c_se <- sqrt(diag(vcovCL(mbL02,cluster=na.omit(dR[,c("start_id",all.vars(mbL02$terms))])$start_id)))
# coeftest(mbL02, vcov.=vcovHC(mbL02, type="HC2"))
# mbL02c_se <- sqrt(diag(vcovHC(mbL02, type="HC2")))
mbL02c_p <- pt(-abs(summary(mbL02)$coefficients[,1]/mbL02c_se), df = mbL02$df.residual)*2
screenreg(list(mbL01,mbL02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbL01c_se,mbL02c_se),
          override.pvalues = list(mbL01c_p,mbL02c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Linear Model (Extended)
#'

## Extended (Linear Model)
mL01 <- lm(limitexpdif ~ word_emperor + word_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dL)
coeftest(mL01, vcov.=vcovCL(mL01,cluster=na.omit(dL[,c("start_id",all.vars(mL01$terms))])$start_id))
mL01c_se <- sqrt(diag(vcovCL(mL01,cluster=na.omit(dL[,c("start_id",all.vars(mL01$terms))])$start_id)))
# coeftest(mL01, vcov.=vcovHC(mL01,type="HC2"))
# mL01c_se <- sqrt(diag(vcovCL(mL01,type="HC2")))
mL01c_p <- pt(-abs(summary(mL01)$coefficients[,1]/mL01c_se), df = mL01$df.residual)*2
mL02 <- lm(limitexpdif ~ word_emperor + word_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dR)
summary(mL02)
coeftest(mL02, vcov.=vcovCL(mL02,cluster=na.omit(dR[,c("start_id",all.vars(mL02$terms))])$start_id))
mL02c_se <- sqrt(diag(vcovCL(mL02,cluster=na.omit(dR[,c("start_id",all.vars(mL02$terms))])$start_id)))
# coeftest(mL02, vcov.=vcovHC(mL02,type="HC2"))
# mL02c_se <- sqrt(diag(vcovCL(mL02,type="HC2")))
mL02c_p <- pt(-abs(summary(mL02)$coefficients[,1]/mL02c_se), df = mL02$df.residual)*2

screenreg(list(mL01,mL02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mL01c_se,mL02c_se),
          override.pvalues = list(mL01c_p,mL02c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Logit (Baseline)
#'

## Baseline (Logit)
mbG01 <- glm(limitexpdec ~ word_emperor + word_emperor_liberal, 
             data = dL, family=binomial("logit"))
coeftest(mbG01, vcov.=vcovCL(mbG01,cluster=na.omit(dL[,c("start_id",all.vars(mbG01$terms))])$start_id))
mbG01c_se <- sqrt(diag(vcovCL(mbG01,cluster=na.omit(dL[,c("start_id",all.vars(mbG01$terms))])$start_id)))
# coeftest(mbG01, vcov.=vcovHC(mbG01, type="HC2"))
# mbG01c_se <- sqrt(diag(vcovHC(mbG01, type="HC2")))
mbG01c_p <- pnorm(-abs(summary(mbG01)$coefficients[,1]/mbG01c_se))*2
mbG02 <- glm(limitexpdec ~ word_emperor + word_emperor_liberal, 
             data = dR, family=binomial("logit"))
coeftest(mbG02, vcov.=vcovCL(mbG02,cluster=na.omit(dR[,c("start_id",all.vars(mbG02$terms))])$start_id))
mbG02c_se <- sqrt(diag(vcovCL(mbG02,cluster=na.omit(dR[,c("start_id",all.vars(mbG02$terms))])$start_id)))
# coeftest(mbG02, vcov.=vcovHC(mbG02, type="HC2"))
# mbG02c_se <- sqrt(diag(vcovHC(mbG02, type="HC2")))
mbG02c_p <- pnorm(-abs(summary(mbG02)$coefficients[,1]/mbG02c_se))*2
screenreg(list(mbG01,mbG02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbG01c_se,mbG02c_se),
          override.pvalues = list(mbG01c_p,mbG02c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Logit (Extended)
#'

## Extended (Logit)
mG01 <- glm(limitexpdec ~ word_emperor + word_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dL, family=binomial("logit"))
coeftest(mG01, vcov.=vcovCL(mG01,cluster=na.omit(dL[,c("start_id",all.vars(mG01$terms))])$start_id))
mG01c_se <- sqrt(diag(vcovCL(mG01,cluster=na.omit(dL[,c("start_id",all.vars(mG01$terms))])$start_id)))
# coeftest(mG01, vcov.=vcovHC(mG01,type="HC2"))
# mG01c_se <- sqrt(diag(vcovCL(mG01,type="HC2")))
mG01c_p <- pnorm(-abs(summary(mG01)$coefficients[,1]/mG01c_se))*2
mG02 <- glm(limitexpdec ~ word_emperor + word_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dR, family=binomial("logit"))
summary(mG02)
coeftest(mG02, vcov.=vcovCL(mG02,cluster=na.omit(dR[,c("start_id",all.vars(mG02$terms))])$start_id))
mG02c_se <- sqrt(diag(vcovCL(mG02,cluster=na.omit(dR[,c("start_id",all.vars(mG02$terms))])$start_id)))
# coeftest(mG02, vcov.=vcovHC(mG02,type="HC2"))
# mG02c_se <- sqrt(diag(vcovCL(mG02,type="HC2")))
mG02c_p <- pnorm(-abs(summary(mG02)$coefficients[,1]/mG02c_se))*2

screenreg(list(mG01,mG02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mG01c_se,mG02c_se),
          override.pvalues = list(mG01c_p,mG02c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Ordinal Logit (Baseline)
#'

## Baseline (Ordinal logit)
mbD01 <- polr(as.ordered(limitexpdifx) ~ word_emperor + word_emperor_liberal, 
            data = dL, Hess = T)
summary(mbD01)
coeftest(mbD01, vcov.=vcovCL(mbD01,cluster=na.omit(dL[,c("start_id",all.vars(mbD01$terms))])$start_id))
mbD01c_se <- sqrt(diag(vcovCL(mbD01,cluster=na.omit(dL[,c("start_id",all.vars(mbD01$terms))])$start_id)))
mbD01c_p <- pnorm(-abs(summary(mbD01)$coefficients[,1]/mbD01c_se))*2
mbD02 <- polr(as.ordered(limitexpdifx) ~ word_emperor + word_emperor_liberal, 
              data = dR, Hess = T)
summary(mbD02)
coeftest(mbD02, vcov.=vcovCL(mbD02,cluster=na.omit(dR[,c("start_id",all.vars(mbD02$terms))])$start_id))
mbD02c_se <- sqrt(diag(vcovCL(mbD02,cluster=na.omit(dR[,c("start_id",all.vars(mbD02$terms))])$start_id)))
mbD02c_p <- pnorm(-abs(summary(mbD02)$coefficients[,1]/mbD02c_se))*2

screenreg(list(mbD01,mbD02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbD01c_se,mbD02c_se),
          override.pvalues = list(mbD01c_p,mbD02c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Ordinal Logit (Extended)
#'

## Extended (Ordinal logit)
mD01 <- polr(as.ordered(limitexpdifx) ~ word_emperor + word_emperor_liberal + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
          data = dL, Hess = T)
summary(mD01)
coeftest(mD01, vcov.=vcovCL(mD01,cluster=na.omit(dL[,c("start_id",all.vars(mD01$terms))])$start_id))
mD01c_se <- sqrt(diag(vcovCL(mD01,cluster=na.omit(dL[,c("start_id",all.vars(mD01$terms))])$start_id)))
mD01c_p <- pnorm(-abs(summary(mD01)$coefficients[,1]/mD01c_se))*2
mD02 <- polr(as.ordered(limitexpdifx) ~ word_emperor + word_emperor_liberal + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
          data = dR, Hess=T)
summary(mD02)
coeftest(mD02, vcov.=vcovCL(mD02,cluster=na.omit(dR[,c("start_id",all.vars(mD02$terms))])$start_id))
mD02c_se <- sqrt(diag(vcovCL(mD02,cluster=na.omit(dR[,c("start_id",all.vars(mD02$terms))])$start_id)))
mD02c_p <- pnorm(-abs(summary(mD02)$coefficients[,1]/mD02c_se))*2

screenreg(list(mD01,mD02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mD01c_se,mD02c_se),
          override.pvalues = list(mD01c_p,mD02c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Ordinal Logit with 3 Category DV (Table G.8)
#'

## 3 Category Model
mC01 <- polr(as.ordered(limitexpdif3) ~ word_emperor + word_emperor_liberal, 
           data = dL, Hess = T)
summary(mC01)
coeftest(mC01, vcov.=vcovCL(mC01,cluster=na.omit(dL[,c("start_id",all.vars(mC01$terms))])$start_id))
mC01c_se <- sqrt(diag(vcovCL(mC01,cluster=na.omit(dL[,c("start_id",all.vars(mC01$terms))])$start_id)))
mC01c_p <- pnorm(-abs(summary(mC01)$coefficients[,1]/mC01c_se))*2
mC02 <- polr(as.ordered(limitexpdif3) ~ word_emperor + word_emperor_liberal,
           data = dR, Hess = T)
summary(mC02)
coeftest(mC02, vcov.=vcovCL(mC02,cluster=na.omit(dR[,c("start_id",all.vars(mC02$terms))])$start_id))
mC02c_se <- sqrt(diag(vcovCL(mC02,cluster=na.omit(dR[,c("start_id",all.vars(mC02$terms))])$start_id)))
mC02c_p <- pnorm(-abs(summary(mC02)$coefficients[,1]/mC02c_se))*2

screenreg(list(mC01,mC02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mC01c_se,mC02c_se),
          override.pvalues = list(mC01c_p,mC02c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mC01,mC02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
       beside = T, digits = 3, single.row = T,
       override.se = list(mC01c_se,mC02c_se),
       override.pvalues = list(mC01c_p,mC02c_p),
       custom.model.names = c("Hate Speech","Biased History"),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "Endorsement treatment effects on the difference in the support for regulating expression in public places (three category outcome, ordinal logit)",
       file = "../out/resout_endorsement_diff_ol_3out_v2.tex", 
       label = "table:resout_endorsement_diff_ol_3out")

#'
#' ### Panel Models
#' 

## Subset data
dLpanel <- subset(dpanel, frame_left==1)
dRpanel <- subset(dpanel, frame_right==1)

#'
#' #### Linear Model (Baseline)
#'

## Baseline (Linear Model)
mbL1 <- lm(limitexp ~ 
            after + Aword_emperor + Aword_emperor_liberal + 
            Bword_emperor + Bword_emperor_liberal, 
          data = dLpanel)
coeftest(mbL1, vcov.=vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id))
mbL1c_se <- sqrt(diag(vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id)))
mbL1c_p <- pt(-abs(summary(mbL1)$coefficients[,1]/mbL1c_se), df = mbL1$df.residual)*2
mbL2 <- lm(limitexp ~ 
            after + Aword_emperor + Aword_emperor_liberal + 
            Bword_emperor + Bword_emperor_liberal, 
          data = dRpanel)
summary(mbL2)
coeftest(mbL2, vcov.=vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id))
mbL2c_se <- sqrt(diag(vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id)))
mbL2c_p <- pt(-abs(summary(mbL2)$coefficients[,1]/mbL2c_se), df = mbL2$df.residual)*2

screenreg(list(mbL1,mbL2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbL1c_se,mbL2c_se),
          override.pvalues = list(mbL1c_p,mbL2c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Linear Model (Extended)
#'

## Extended (Linear Model)
mL1 <- lm(limitexp ~ 
            after + Aword_emperor + Aword_emperor_liberal + 
            Bword_emperor + Bword_emperor_liberal + 
            fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
          data = dLpanel)
summary(mL1)
coeftest(mL1, vcov.=vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id))
mL1c_se <- sqrt(diag(vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id)))
mL1c_p <- pt(-abs(summary(mL1)$coefficients[,1]/mL1c_se), df = mL1$df.residual)*2
mL2 <- lm(limitexp ~ 
            after + Aword_emperor + Aword_emperor_liberal + 
            Bword_emperor + Bword_emperor_liberal + 
            fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
          data = dRpanel)
summary(mL2)
coeftest(mL2, vcov.=vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id))
mL2c_se <- sqrt(diag(vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id)))
mL2c_p <- pt(-abs(summary(mL2)$coefficients[,1]/mL2c_se), df = mL2$df.residual)*2

screenreg(list(mL1,mL2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mL1c_se,mL2c_se),
          override.pvalues = list(mL1c_p,mL2c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Logit (Baseline)
#'

## Baseline (Logit)
mbG1 <- glm(limitexp > 1 ~ 
             after + Aword_emperor + Aword_emperor_liberal + 
             Bword_emperor + Bword_emperor_liberal, 
           data = dLpanel, family=binomial("logit"))
coeftest(mbG1, vcov.=vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id))
mbG1c_se <- sqrt(diag(vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id)))
mbG1c_p <- pnorm(-abs(summary(mbG1)$coefficients[,1]/mbG1c_se))*2
mbG2 <- glm(limitexp > 1 ~ 
             after + Aword_emperor + Aword_emperor_liberal + 
             Bword_emperor + Bword_emperor_liberal, 
           data = dRpanel, family=binomial("logit"))
summary(mbG2)
coeftest(mbG2, vcov.=vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id))
mbG2c_se <- sqrt(diag(vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id)))
mbG2c_p <- pnorm(-abs(summary(mbG2)$coefficients[,1]/mbG2c_se))*2

screenreg(list(mbG1,mbG2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbG1c_se,mbG2c_se),
          override.pvalues = list(mbG1c_p,mbG2c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Logit (Extended)
#'

## Extended (Logit)
mG1 <- glm(limitexp > 1 ~ 
            after + Aword_emperor + Aword_emperor_liberal + 
            Bword_emperor + Bword_emperor_liberal + 
            fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
          data = dLpanel, family=binomial("logit"))
summary(mG1)
coeftest(mG1, vcov.=vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id))
mG1c_se <- sqrt(diag(vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id)))
mG1c_p <- pnorm(-abs(summary(mG1)$coefficients[,1]/mG1c_se))*2
mG2 <- glm(limitexp > 1 ~ 
            after + Aword_emperor + Aword_emperor_liberal + 
            Bword_emperor + Bword_emperor_liberal + 
            fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
          data = dRpanel, family=binomial("logit"))
summary(mG2)
coeftest(mG2, vcov.=vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id))
mG2c_se <- sqrt(diag(vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id)))
mG2c_p <- pnorm(-abs(summary(mG2)$coefficients[,1]/mG2c_se))*2

screenreg(list(mG1,mG2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mG1c_se,mG2c_se),
          override.pvalues = list(mG1c_p,mG2c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Ordinal Logit (Baseline)
#'

## Baseline (Ordinal logit)
mbD1 <- polr(as.ordered(limitexp) ~ 
               after + Aword_emperor + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal, 
             data = dLpanel, Hess = T)
summary(mbD1)
coeftest(mbD1, vcov.=vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id))
mbD1c_se <- sqrt(diag(vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id)))
mbD1c_p <- pnorm(-abs(summary(mbD1)$coefficients[,1]/mbD1c_se))*2

mbD2 <- polr(as.ordered(limitexp) ~ 
               after + Aword_emperor + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal, 
             data = dRpanel, Hess = T)
summary(mbD2)
coeftest(mbD2, vcov.=vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id))
mbD2c_se <- sqrt(diag(vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id)))
mbD2c_p <- pnorm(-abs(summary(mbD2)$coefficients[,1]/mbD2c_se))*2

screenreg(list(mbD1,mbD2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mbD1c_se,mbD2c_se),
          override.pvalues = list(mbD1c_p,mbD2c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' #### Ordinal Logit (Extended)
#'

## Extended (Ordinal logit)
mD1 <- polr(as.ordered(limitexp) ~ 
             after + Aword_emperor + Aword_emperor_liberal + 
             Bword_emperor + Bword_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dLpanel, Hess = T)
summary(mD1)
coeftest(mD1, vcov.=vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id))
mD1c_se <- sqrt(diag(vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id)))
mD1c_p <- pnorm(-abs(summary(mD1)$coefficients[,1]/mD1c_se))*2

mD2 <- polr(as.ordered(limitexp) ~ 
             after + Aword_emperor + Aword_emperor_liberal + 
             Bword_emperor + Bword_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dRpanel, Hess = T)
summary(mD2)
coeftest(mD2, vcov.=vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id))
mD2c_se <- sqrt(diag(vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id)))
mD2c_p <- pnorm(-abs(summary(mD2)$coefficients[,1]/mD2c_se))*2

screenreg(list(mD1,mD2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = T,
          override.se = list(mD1c_se,mD2c_se),
          override.pvalues = list(mD1c_p,mD2c_p),
          custom.model.names = c("Hate Speech","Biased History"),
          custom.coef.map = vn)

#'
#' ### Combined Tables
#'
#' #### Linear Model (Baseline) (Table G.2)
#'

## Baseline (Linear Regression)
screenreg(list(mbL1,mbL01,mbL2,mbL02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = F,
          override.se = list(mbL1c_se,mbL01c_se,mbL2c_se,mbL02c_se),
          override.pvalues = list(mbL1c_p,mbL01c_p,mbL2c_p,mbL02c_p),
          custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
          custom.model.names = rep(c("Panel","Difference"),2),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mbL1,mbL01,mbL2,mbL02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
       beside = T, digits = 3, single.row = F,
       override.se = list(mbL1c_se,mbL01c_se,mbL2c_se,mbL02c_se),
       override.pvalues = list(mbL1c_p,mbL01c_p,mbL2c_p,mbL02c_p),
       custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
       custom.model.names = rep(c("Panel","Difference"),2),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "Endorsement treatment effects on the support for regulating expression in public places (OLS)",
       file = "../out/resout_endorsement_lm_base_v2.tex", 
       label = "table:resout_endorsement_lm_base")

#'
#' #### Linear Model (Extended) (Table G.3)
#'

## Extended (Linear Regression)
screenreg(list(mL1,mL01,mL2,mL02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = F,
          override.se = list(mL1c_se,mL01c_se,mL2c_se,mL02c_se),
          override.pvalues = list(mL1c_p,mL01c_p,mL2c_p,mL02c_p),
          custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
          custom.model.names = rep(c("Panel","Difference"),2),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mL1,mL01,mL2,mL02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
       beside = T, digits = 3, single.row = F,
       override.se = list(mL1c_se,mL01c_se,mL2c_se,mL02c_se),
       override.pvalues = list(mL1c_p,mL01c_p,mL2c_p,mL02c_p),
       custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
       custom.model.names = rep(c("Panel","Difference"),2),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "Endorsement treatment effects on the support for regulating expression in public places (OLS, extended models)",
       file = "../out/resout_endorsement_lm_ext_v2.tex", 
       label = "table:resout_endorsement_lm_ext")
#+

#'
#' ### Logit (Baseline) (Table G.4)
#'

## Baseline (Logit)
screenreg(list(mbG1,mbG01,mbG2,mbG02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = F,
          override.se = list(mbG1c_se,mbG01c_se,mbG2c_se,mbG02c_se),
          override.pvalues = list(mbG1c_p,mbG01c_p,mbG2c_p,mbG02c_p),
          custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
          custom.model.names = rep(c("Panel","Difference"),2),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mbG1,mbG01,mbG2,mbG02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
       beside = T, digits = 3, single.row = F,
       override.se = list(mbG1c_se,mbG01c_se,mbG2c_se,mbG02c_se),
       override.pvalues = list(mbG1c_p,mbG01c_p,mbG2c_p,mbG02c_p),
       custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
       custom.model.names = rep(c("Panel","Difference"),2),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "Endorsement treatment effects on the support for regulating expression in public places (Logit)",
       file = "../out/resout_endorsement_logit_base_v2.tex", 
       label = "table:resout_endorsement_logit_base")

#'
#' #### Logit (Extended) (Table G.5)
#'

## Extended (Logit)
screenreg(list(mG1,mG01,mG2,mG02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = F,
          override.se = list(mG1c_se,mG01c_se,mG2c_se,mG02c_se),
          override.pvalues = list(mG1c_p,mG01c_p,mG2c_p,mG02c_p),
          custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
          custom.model.names = rep(c("Panel","Difference"),2),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mG1,mG01,mG2,mG02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
       beside = T, digits = 3, single.row = F,
       override.se = list(mG1c_se,mG01c_se,mG2c_se,mG02c_se),
       override.pvalues = list(mG1c_p,mG01c_p,mG2c_p,mG02c_p),
       custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
       custom.model.names = rep(c("Panel","Difference"),2),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "Endorsement treatment effects on the support for regulating expression in public places (Logit, extended models)",
       file = "../out/resout_endorsement_logit_ext_v2.tex", 
       label = "table:resout_endorsement_logit_ext")
#+ 

#'
#' ### Ordinal Logit (Baseline) (Table G.6)
#'

## Baseline (Ordinal logit)
screenreg(list(mbD1,mbD01,mbD2,mbD02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = F,
          override.se = list(mbD1c_se,mbD01c_se,mbD2c_se,mbD02c_se),
          override.pvalues = list(mbD1c_p,mbD01c_p,mbD2c_p,mbD02c_p),
          custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
          custom.model.names = rep(c("Panel","Difference"),2),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mbD1,mbD01,mbD2,mbD02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
       beside = T, digits = 3, single.row = F,
       override.se = list(mbD1c_se,mbD01c_se,mbD2c_se,mbD02c_se),
       override.pvalues = list(mbD1c_p,mbD01c_p,mbD2c_p,mbD02c_p),
       custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
       custom.model.names = rep(c("Panel","Difference"),2),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "Endorsement treatment effects on the support for regulating expression in public places (ordinal logit)",
       file = "../out/resout_endorsement_ol_base_v2.tex",
       label = "table:resout_endorsement_ol_base")

#'
#' #### Ordinal Logit (Extended) (Table G.7)
#'

## Extended (Ordinal logit)
screenreg(list(mD1,mD01,mD2,mD02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
          beside = T, digits = 3, single.row = F,
          override.se = list(mD1c_se,mD01c_se,mD2c_se,mD02c_se),
          override.pvalues = list(mD1c_p,mD01c_p,mD2c_p,mD02c_p),
          custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
          custom.model.names = rep(c("Panel","Difference"),2),
          custom.coef.map = vn)

#+ eval=FALSE
texreg(list(mD1,mD01,mD2,mD02), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
       beside = T, digits = 3, single.row = F,
       override.se = list(mD1c_se,mD01c_se,mD2c_se,mD02c_se),
       override.pvalues = list(mD1c_p,mD01c_p,mD2c_p,mD02c_p),
       custom.header = list("Hate Speech" = 1:2, "Biased History" = 3:4),
       custom.model.names = rep(c("Panel","Difference"),2),
       custom.coef.map = vn,
       custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       caption = "Endorsement treatment effects on the support for regulating expression in public places (ordinal logit, extended model)",
       file = "../out/resout_endorsement_ol_ext_v2.tex",
       label = "table:resout_endorsement_ol_ext")

#'
#' ### Plotting Effects
#' 
#' #### Linear Model (Baseline) (Figure 2)
#' 

## Linear Regression (Baseline) ##

tmp1 <- lm(limitexp ~ 
              after + Aword_expert + Aword_emperor_liberal + 
              Bword_expert + Bword_emperor_liberal, 
            data = dLpanel)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp ~ 
               after + Aword_emperor + Aword_expert + 
               Bword_emperor + Bword_expert, 
             data = dLpanel)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt1 <- 
  rbind(c(mbL1$coefficients[2],mbL1c_se[2],
          coefci(mbL1, vcov.=vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id), level=0.95)[2,],
          coefci(mbL1, vcov.=vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt1 <- as.data.frame(tmpdt1)
colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt1$frame <- "Hate speech"
tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")

tmp1 <- lm(limitexp ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_expert + Bword_emperor_liberal, 
           data = dRpanel)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_expert, 
           data = dRpanel)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt2 <- 
  rbind(c(mbL2$coefficients[2],mbL2c_se[2],
          coefci(mbL2, vcov.=vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id), level=0.95)[2,],
          coefci(mbL2, vcov.=vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt2 <- as.data.frame(tmpdt2)
colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt2$frame <- "Biased history"
tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")

## Combine Data
mbL1dt <- rbind(tmpdt1,tmpdt2)
mbL1dt$frame <- factor(mbL1dt$frame, levels = c("Hate speech","Biased history"))
mbL1dt$endorsement <- factor(mbL1dt$endorsement, levels = unique(mbL1dt$endorsement))

write.csv(mbL1dt, row.names = F, file = "../out/effect_endorsement_lm_base_v2.csv")

## Plot of Endorsement Treatment Effects

require(ggplot2)

p <- ggplot(mbL1dt, aes(x=endorsement,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~frame) + 
  labs(x="Endorser", y="Negative endorsement treatment effect\nOLS coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_endorsement_lm_base_v2.pdf", width=6, height=4)
ggsave("../out/effect_endorsement_lm_base_v2.png", width=6, height=4)

#' 
#' #### Linear Model (Extended) (Figure G.1)
#' 

## Linear Regression (Extended) ##

tmp1 <- lm(limitexp ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_expert + Bword_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dLpanel)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_expert + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dLpanel)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt1 <- 
  rbind(c(mL1$coefficients[2],mL1c_se[2],
          coefci(mL1, vcov.=vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id), level=0.95)[2,],
          coefci(mL1, vcov.=vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt1 <- as.data.frame(tmpdt1)
colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt1$frame <- "Hate speech"
tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")

tmp1 <- lm(limitexp ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_expert + Bword_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dRpanel)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
tmp2 <- lm(limitexp ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_expert + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dRpanel)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2

tmpdt2 <- 
  rbind(c(mL2$coefficients[2],mL2c_se[2],
          coefci(mL2, vcov.=vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id), level=0.95)[2,],
          coefci(mL2, vcov.=vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt2 <- as.data.frame(tmpdt2)
colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt2$frame <- "Biased history"
tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")

## Combine Data
mL1dt <- rbind(tmpdt1,tmpdt2)
mL1dt$frame <- factor(mL1dt$frame, levels = c("Hate speech","Biased history"))
mL1dt$endorsement <- factor(mL1dt$endorsement, levels = unique(mL1dt$endorsement))

write.csv(mL1dt, row.names = F, file = "../out/effect_endorsement_lm_ext_v2.csv")

## Plot of Endorsement Treatment Effects

require(ggplot2)

p <- ggplot(mL1dt, aes(x=endorsement,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~frame) + 
  labs(x="Endorser", y="Negative endorsement treatment effect\nOLS coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_endorsement_lm_ext_v2.pdf", width=6, height=4)
ggsave("../out/effect_endorsement_lm_ext_v2.png", width=6, height=4)

#' 
#' #### Logit (Baseline) (Figure G.2)
#' 

## Logit (Baseline) ##

tmp1 <- glm(limitexp > 1 ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_expert + Bword_emperor_liberal, 
           data = dLpanel, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp > 1 ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_expert, 
           data = dLpanel, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt1 <- 
  rbind(c(mbG1$coefficients[2],mbG1c_se[2],
          coefci(mbG1, vcov.=vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id), level=0.95)[2,],
          coefci(mbG1, vcov.=vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt1 <- as.data.frame(tmpdt1)
colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt1$frame <- "Hate speech"
tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")

tmp1 <- glm(limitexp > 1 ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_expert + Bword_emperor_liberal, 
           data = dRpanel, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp > 1 ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_expert, 
           data = dRpanel, family=binomial("logit"))
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt2 <- 
  rbind(c(mbG2$coefficients[2],mbG2c_se[2],
          coefci(mbG2, vcov.=vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id), level=0.95)[2,],
          coefci(mbG2, vcov.=vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt2 <- as.data.frame(tmpdt2)
colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt2$frame <- "Biased history"
tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")

## Combine Data
mbG1dt <- rbind(tmpdt1,tmpdt2)
mbG1dt$frame <- factor(mbG1dt$frame, levels = c("Hate speech","Biased history"))
mbG1dt$endorsement <- factor(mbG1dt$endorsement, levels = unique(mbG1dt$endorsement))

write.csv(mbG1dt, row.names = F, file = "../out/effect_endorsement_logit_base_v2.csv")

## Plot of Endorsement Treatment Effects

require(ggplot2)

p <- ggplot(mbG1dt, aes(x=endorsement,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~frame) + 
  labs(x="Endorser", y="Negative endorsement treatment effect\nLogit coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_endorsement_logit_base_v2.pdf", width=6, height=4)
ggsave("../out/effect_endorsement_logit_base_v2.png", width=6, height=4)

#' 
#' #### Logit (Extended) (Figure G.3)
#' 

## Logit (Extended) ##

tmp1 <- glm(limitexp > 1 ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_expert + Bword_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dLpanel, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_expert + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dLpanel)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt1 <- 
  rbind(c(mG1$coefficients[2],mG1c_se[2],
          coefci(mG1, vcov.=vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id), level=0.95)[2,],
          coefci(mG1, vcov.=vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt1 <- as.data.frame(tmpdt1)
colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt1$frame <- "Hate speech"
tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")

tmp1 <- glm(limitexp > 1 ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_expert + Bword_emperor_liberal + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dRpanel, family=binomial("logit"))
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- glm(limitexp ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_expert + 
             fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
           data = dRpanel)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt2 <- 
  rbind(c(mG2$coefficients[2],mG2c_se[2],
          coefci(mG2, vcov.=vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id), level=0.95)[2,],
          coefci(mG2, vcov.=vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id), level=0.90)[2,]),
        c(tmp1$coefficients[2],tmp1c_se[2],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
        c(tmp2$coefficients[2],tmp2c_se[2],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))

tmpdt2 <- as.data.frame(tmpdt2)
colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt2$frame <- "Biased history"
tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")

## Combine Data
mG1dt <- rbind(tmpdt1,tmpdt2)
mG1dt$frame <- factor(mG1dt$frame, levels = c("Hate speech","Biased history"))
mG1dt$endorsement <- factor(mG1dt$endorsement, levels = unique(mG1dt$endorsement))

write.csv(mG1dt, row.names = F, file = "../out/effect_endorsement_logit_ext_v2.csv")

## Plot of Endorsement Treatment Effects

require(ggplot2)

p <- ggplot(mG1dt, aes(x=endorsement,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~frame) + 
  labs(x="Endorser", y="Negative endorsement treatment effect\nLogit coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_endorsement_logit_ext_v2.pdf", width=6, height=4)
ggsave("../out/effect_endorsement_logit_ext_v2.png", width=6, height=4)

#' 
#' #### Ordinal Logit (Baseline) (Figure G.4)
#' 

## Ordinal Logit (Baseline) ##

tmp1 <- polr(as.ordered(limitexp) ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_emperor + Bword_emperor_liberal, 
           data = dLpanel, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp) ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_emperor_liberal, 
           data = dLpanel, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt1 <- 
  rbind(c(mbD1$coefficients[1],mbD1c_se[1],
          coefci(mbD1, vcov.=vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id), level=0.95)[1,],
          coefci(mbD1, vcov.=vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt1 <- as.data.frame(tmpdt1)
colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt1$frame <- "Hate speech"
tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")

tmp1 <- polr(as.ordered(limitexp) ~ 
             after + Aword_expert + Aword_emperor_liberal + 
             Bword_emperor + Bword_emperor_liberal, 
           data = dRpanel, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp) ~ 
             after + Aword_emperor + Aword_expert + 
             Bword_emperor + Bword_emperor_liberal, 
           data = dRpanel, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt2 <- 
  rbind(c(mbD2$coefficients[1],mbD2c_se[1],
          coefci(mbD2, vcov.=vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id), level=0.95)[1,],
          coefci(mbD2, vcov.=vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt2 <- as.data.frame(tmpdt2)
colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt2$frame <- "Biased history"
tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")

## Combine Data
mbD1dt <- rbind(tmpdt1,tmpdt2)
mbD1dt$frame <- factor(mbD1dt$frame, levels = c("Hate speech","Biased history"))
mbD1dt$endorsement <- factor(mbD1dt$endorsement, levels = unique(mbD1dt$endorsement))

write.csv(mbD1dt, row.names = F, file = "../out/effect_endorsement_ol_base_v2.csv")

## Plot of Endorsement Treatment Effects

require(ggplot2)

p <- ggplot(mbD1dt, aes(x=endorsement,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~frame) + 
  labs(x="Endorser", y="Negative endorsement treatment effect\nordinal logit coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_endorsement_ol_base_v2.pdf", width=6, height=4)
ggsave("../out/effect_endorsement_ol_base_v2.png", width=6, height=4)

#' 
#' #### Ordinal Logit (Extended) (Figure G.5)
#' 

## Ordinal Logit (Extended) ##

tmp1 <- polr(as.ordered(limitexp) ~ 
               after + Aword_expert + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dLpanel, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp) ~ 
               after + Aword_emperor + Aword_expert + 
               Bword_emperor + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dLpanel, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt1 <- 
  rbind(c(mD1$coefficients[1],mD1c_se[1],
          coefci(mD1, vcov.=vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id), level=0.95)[1,],
          coefci(mD1, vcov.=vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt1 <- as.data.frame(tmpdt1)
colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt1$frame <- "Hate speech"
tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")

tmp1 <- polr(as.ordered(limitexp) ~ 
               after + Aword_expert + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dRpanel, Hess = TRUE)
summary(tmp1)
coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
tmp2 <- polr(as.ordered(limitexp) ~ 
               after + Aword_emperor + Aword_expert + 
               Bword_emperor + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dRpanel, Hess = TRUE)
summary(tmp2)
coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2

tmpdt2 <- 
  rbind(c(mD2$coefficients[1],mD2c_se[1],
          coefci(mD2, vcov.=vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id), level=0.95)[1,],
          coefci(mD2, vcov.=vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id), level=0.90)[1,]),
        c(tmp1$coefficients[1],tmp1c_se[1],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
          coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
        c(tmp2$coefficients[1],tmp2c_se[1],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
          coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))

tmpdt2 <- as.data.frame(tmpdt2)
colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
tmpdt2$frame <- "Biased history"
tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")

## Combine Data
mD1dt <- rbind(tmpdt1,tmpdt2)
mD1dt$frame <- factor(mD1dt$frame, levels = c("Hate speech","Biased history"))
mD1dt$endorsement <- factor(mD1dt$endorsement, levels = unique(mD1dt$endorsement))

write.csv(mD1dt, row.names = F, file = "../out/effect_endorsement_ol_ext_v2.csv")

## Plot of Endorsement Treatment Effects

require(ggplot2)

p <- ggplot(mD1dt, aes(x=endorsement,y=cf)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.5) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2) + 
  geom_point(size=2) + 
  facet_grid(.~frame) + 
  labs(x="Endorser", y="Negative endorsement treatment effect\nordinal logit coefficient") + 
  theme_bw()

#+ fig.width=6, fig.height=4
p

#+ eval=FALSE
ggsave("../out/effect_endorsement_ol_ext_v2.pdf", width=6, height=4)
ggsave("../out/effect_endorsement_ol_ext_v2.png", width=6, height=4)


#'
#' ## Moderation by Ideologues vs. Moderates
#'

tmpfunc <- function(dLpanel, dRpanel, postfixset) {
  
  ## Baseline (Linear Model)
  mbL1 <- lm(limitexp ~ 
               after + Aword_emperor + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal, 
             data = dLpanel)
  coeftest(mbL1, vcov.=vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id))
  mbL1c_se <- sqrt(diag(vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id)))
  mbL1c_p <- pt(-abs(summary(mbL1)$coefficients[,1]/mbL1c_se), df = mbL1$df.residual)*2
  mbL2 <- lm(limitexp ~ 
               after + Aword_emperor + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal, 
             data = dRpanel)
  summary(mbL2)
  coeftest(mbL2, vcov.=vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id))
  mbL2c_se <- sqrt(diag(vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id)))
  mbL2c_p <- pt(-abs(summary(mbL2)$coefficients[,1]/mbL2c_se), df = mbL2$df.residual)*2
  
  ## Extended (Linear Model)
  mL1 <- lm(limitexp ~ 
              after + Aword_emperor + Aword_emperor_liberal + 
              Bword_emperor + Bword_emperor_liberal + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
            data = dLpanel)
  summary(mL1)
  coeftest(mL1, vcov.=vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id))
  mL1c_se <- sqrt(diag(vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id)))
  mL1c_p <- pt(-abs(summary(mL1)$coefficients[,1]/mL1c_se), df = mL1$df.residual)*2
  mL2 <- lm(limitexp ~ 
              after + Aword_emperor + Aword_emperor_liberal + 
              Bword_emperor + Bword_emperor_liberal + 
              fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
            data = dRpanel)
  summary(mL2)
  coeftest(mL2, vcov.=vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id))
  mL2c_se <- sqrt(diag(vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id)))
  mL2c_p <- pt(-abs(summary(mL2)$coefficients[,1]/mL2c_se), df = mL2$df.residual)*2
  
  ## Baseline (logit)
  mbG1 <- glm(limitexp > 1 ~ 
                after + Aword_emperor + Aword_emperor_liberal + 
                Bword_emperor + Bword_emperor_liberal, 
              data = dLpanel, family=binomial("logit"))
  coeftest(mbG1, vcov.=vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id))
  mbG1c_se <- sqrt(diag(vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id)))
  mbG1c_p <- pnorm(-abs(summary(mbG1)$coefficients[,1]/mbG1c_se))*2
  mbG2 <- glm(limitexp > 1 ~ 
                after + Aword_emperor + Aword_emperor_liberal + 
                Bword_emperor + Bword_emperor_liberal, 
              data = dRpanel, family=binomial("logit"))
  summary(mbG2)
  coeftest(mbG2, vcov.=vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id))
  mbG2c_se <- sqrt(diag(vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id)))
  mbG2c_p <- pnorm(-abs(summary(mbG2)$coefficients[,1]/mbG2c_se))*2
  
  ## Extended (logit)
  mG1 <- glm(limitexp > 1 ~ 
               after + Aword_emperor + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dLpanel, family=binomial("logit"))
  summary(mG1)
  coeftest(mG1, vcov.=vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id))
  mG1c_se <- sqrt(diag(vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id)))
  mG1c_p <- pnorm(-abs(summary(mG1)$coefficients[,1]/mG1c_se))*2
  mG2 <- glm(limitexp > 1 ~ 
               after + Aword_emperor + Aword_emperor_liberal + 
               Bword_emperor + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dRpanel, family=binomial("logit"))
  summary(mG2)
  coeftest(mG2, vcov.=vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id))
  mG2c_se <- sqrt(diag(vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id)))
  mG2c_p <- pnorm(-abs(summary(mG2)$coefficients[,1]/mG2c_se))*2
  
  ## Baseline (ordinal logit)
  mbD1 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_emperor + Aword_emperor_liberal + 
                 Bword_emperor + Bword_emperor_liberal, 
               data = dLpanel, Hess = T)
  summary(mbD1)
  coeftest(mbD1, vcov.=vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id))
  mbD1c_se <- sqrt(diag(vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id)))
  mbD1c_p <- pnorm(-abs(summary(mbD1)$coefficients[,1]/mbD1c_se))*2
  
  mbD2 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_emperor + Aword_emperor_liberal + 
                 Bword_emperor + Bword_emperor_liberal, 
               data = dRpanel, Hess = T)
  summary(mbD2)
  coeftest(mbD2, vcov.=vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id))
  mbD2c_se <- sqrt(diag(vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id)))
  mbD2c_p <- pnorm(-abs(summary(mbD2)$coefficients[,1]/mbD2c_se))*2
  
  ## Extended (ordinal logit)
  mD1 <- polr(as.ordered(limitexp) ~ 
                after + Aword_emperor + Aword_emperor_liberal + 
                Bword_emperor + Bword_emperor_liberal + 
                fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
              data = dLpanel, Hess = T)
  summary(mD1)
  coeftest(mD1, vcov.=vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id))
  mD1c_se <- sqrt(diag(vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id)))
  mD1c_p <- pnorm(-abs(summary(mD1)$coefficients[,1]/mD1c_se))*2
  
  mD2 <- polr(as.ordered(limitexp) ~ 
                after + Aword_emperor + Aword_emperor_liberal + 
                Bword_emperor + Bword_emperor_liberal + 
                fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
              data = dRpanel, Hess = T)
  summary(mD2)
  coeftest(mD2, vcov.=vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id))
  mD2c_se <- sqrt(diag(vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id)))
  mD2c_p <- pnorm(-abs(summary(mD2)$coefficients[,1]/mD2c_se))*2
  
  ## Combined Tables

  ## Baseline (Linear Regression)
  texreg(list(mbL1,mbL2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
         beside = T, digits = 3, single.row = T,
         override.se = list(mbL1c_se,mbL2c_se),
         override.pvalues = list(mbL1c_p,mbL2c_p),
         custom.model.names = c("Hate Speech", "Biased History"),
         custom.coef.map = vn,
         custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
         use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
         caption = "Endorsement treatment effects on the support for regulating expression in public places (OLS)",
         file = paste0("../out/resout_endorsement_lm_base_v2","_",postfixset,".tex"), 
         label = paste0("table:resout_endorsement_lm_base","_",postfixset))

  ## Extended (Linear Regression)
  texreg(list(mL1,mL2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
         beside = T, digits = 3, single.row = T,
         override.se = list(mL1c_se,mL2c_se),
         override.pvalues = list(mL1c_p,mL2c_p),
         custom.model.names = c("Hate Speech", "Biased History"),
         custom.coef.map = vn,
         custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
         use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
         caption = "Endorsement treatment effects on the support for regulating expression in public places (OLS, extended models)",
         file = paste0("../out/resout_endorsement_lm_ext_v2","_",postfixset,".tex"), 
         label = paste0("table:resout_endorsement_lm_ext","_",postfixset))
  
  ## Baseline (logit)
  texreg(list(mbG1,mbG2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
         beside = T, digits = 3, single.row = F,
         override.se = list(mbG1c_se,mbG2c_se),
         override.pvalues = list(mbG1c_p,mbG2c_p),
         custom.model.names = c("Hate Speech", "Biased History"),
         custom.coef.map = vn,
         custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
         use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
         caption = "Endorsement treatment effects on the support for regulating expression in public places (logit)",
         file = paste0("../out/resout_endorsement_logit_base_v2","_",postfixset,".tex"), 
         label = paste0("table:resout_endorsement_logit_base","_",postfixset))
  
  ## Extended (logit)
  texreg(list(mG1,mG2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
         beside = T, digits = 3, single.row = F,
         override.se = list(mG1c_se,mG2c_se),
         override.pvalues = list(mG1c_p,mG2c_p),
         custom.model.names = c("Hate Speech", "Biased History"),
         custom.coef.map = vn,
         custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
         use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
         caption = "Endorsement treatment effects on the support for regulating expression in public places (logit, extended models)",
         file = paste0("../out/resout_endorsement_logit_ext_v2","_",postfixset,".tex"), 
         label = paste0("table:resout_endorsement_logit_ext","_",postfixset))

  ## Baseline (ordinal logit)
  texreg(list(mbD1,mbD2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
         beside = T, digits = 3, single.row = F,
         override.se = list(mbD1c_se,mbD2c_se),
         override.pvalues = list(mbD1c_p,mbD2c_p),
         custom.model.names = c("Hate Speech", "Biased History"),
         custom.coef.map = vn,
         custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
         use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
         caption = "Endorsement treatment effects on the support for regulating expression in public places (ordinal logit)",
         file = paste0("../out/resout_endorsement_ol_base_v2","_",postfixset,".tex"), 
         label = paste0("table:resout_endorsement_ol_base","_",postfixset))

  ## Extended (ordinal logit)
  texreg(list(mD1,mD2), stars = c(0.1,0.05,0.01,0.001), symbol = "+", 
         beside = T, digits = 3, single.row = F,
         override.se = list(mD1c_se,mD2c_se),
         override.pvalues = list(mD1c_p,mD2c_p),
         custom.model.names = c("Hate Speech", "Biased History"),
         custom.coef.map = vn,
         custom.note = "%stars. Robust standard errors clustered by respondent ID in parentheses.",
         use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
         caption = "Endorsement treatment effects on the support for regulating expression in public places (ordinal logit, extended model)",
         file = paste0("../out/resout_endorsement_ol_ext_v2","_",postfixset,".tex"), 
         label = paste0("table:resout_endorsement_ol_ext","_",postfixset))
  
  ## Plotting Endorsement Treatment Effects

  ## Linear Regression (Baseline) ##
  
  tmp1 <- lm(limitexp ~ 
               after + Aword_expert + Aword_emperor_liberal + 
               Bword_expert + Bword_emperor_liberal, 
             data = dLpanel)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
  tmp2 <- lm(limitexp ~ 
               after + Aword_emperor + Aword_expert + 
               Bword_emperor + Bword_expert, 
             data = dLpanel)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2
  
  tmpdt1 <- 
    rbind(c(mbL1$coefficients[2],mbL1c_se[2],
            coefci(mbL1, vcov.=vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id), level=0.95)[2,],
            coefci(mbL1, vcov.=vcovCL(mbL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbL1$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt1 <- as.data.frame(tmpdt1)
  colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt1$frame <- "Hate speech"
  tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  tmp1 <- lm(limitexp ~ 
               after + Aword_expert + Aword_emperor_liberal + 
               Bword_expert + Bword_emperor_liberal, 
             data = dRpanel)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
  tmp2 <- lm(limitexp ~ 
               after + Aword_emperor + Aword_expert + 
               Bword_emperor + Bword_expert, 
             data = dRpanel)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2
  
  tmpdt2 <- 
    rbind(c(mbL2$coefficients[2],mbL2c_se[2],
            coefci(mbL2, vcov.=vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id), level=0.95)[2,],
            coefci(mbL2, vcov.=vcovCL(mbL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbL2$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt2 <- as.data.frame(tmpdt2)
  colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt2$frame <- "Biased history"
  tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  ## Combine Data
  mbL1dt <- rbind(tmpdt1,tmpdt2)
  mbL1dt$frame <- factor(mbL1dt$frame, levels = c("Hate speech","Biased history"))
  mbL1dt$endorsement <- factor(mbL1dt$endorsement, levels = unique(mbL1dt$endorsement))
  
  write.csv(mbL1dt, row.names = F, 
            file = paste0("../out/effect_endorsement_lm_base_v2","_",postfixset,".csv"))

  ## Linear Regression (Extended) ##
  
  tmp1 <- lm(limitexp ~ 
               after + Aword_expert + Aword_emperor_liberal + 
               Bword_expert + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dLpanel)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
  tmp2 <- lm(limitexp ~ 
               after + Aword_emperor + Aword_expert + 
               Bword_emperor + Bword_expert + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dLpanel)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2
  
  tmpdt1 <- 
    rbind(c(mL1$coefficients[2],mL1c_se[2],
            coefci(mL1, vcov.=vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id), level=0.95)[2,],
            coefci(mL1, vcov.=vcovCL(mL1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mL1$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt1 <- as.data.frame(tmpdt1)
  colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt1$frame <- "Hate speech"
  tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  tmp1 <- lm(limitexp ~ 
               after + Aword_expert + Aword_emperor_liberal + 
               Bword_expert + Bword_emperor_liberal + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dRpanel)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pt(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se), df = tmp1$df.residual)*2
  tmp2 <- lm(limitexp ~ 
               after + Aword_emperor + Aword_expert + 
               Bword_emperor + Bword_expert + 
               fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
             data = dRpanel)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pt(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se), df = tmp2$df.residual)*2
  
  tmpdt2 <- 
    rbind(c(mL2$coefficients[2],mL2c_se[2],
            coefci(mL2, vcov.=vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id), level=0.95)[2,],
            coefci(mL2, vcov.=vcovCL(mL2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mL2$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt2 <- as.data.frame(tmpdt2)
  colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt2$frame <- "Biased history"
  tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  ## Combine Data
  mL1dt <- rbind(tmpdt1,tmpdt2)
  mL1dt$frame <- factor(mL1dt$frame, levels = c("Hate speech","Biased history"))
  mL1dt$endorsement <- factor(mL1dt$endorsement, levels = unique(mL1dt$endorsement))
  
  write.csv(mL1dt, row.names = F, 
            file = paste0("../out/effect_endorsement_lm_ext_v2","_",postfixset,".csv"))
  
  ## logit (Baseline) ##
  
  tmp1 <- glm(limitexp > 1 ~ 
                after + Aword_expert + Aword_emperor_liberal + 
                Bword_expert + Bword_emperor_liberal, 
              data = dLpanel, family=binomial("logit"))
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- glm(limitexp > 1 ~ 
                after + Aword_emperor + Aword_expert + 
                Bword_emperor + Bword_expert, 
              data = dLpanel, family=binomial("logit"))
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt1 <- 
    rbind(c(mbG1$coefficients[2],mbG1c_se[2],
            coefci(mbG1, vcov.=vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id), level=0.95)[2,],
            coefci(mbG1, vcov.=vcovCL(mbG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbG1$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt1 <- as.data.frame(tmpdt1)
  colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt1$frame <- "Hate speech"
  tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  tmp1 <- glm(limitexp > 1 ~ 
                after + Aword_expert + Aword_emperor_liberal + 
                Bword_expert + Bword_emperor_liberal, 
              data = dRpanel, family=binomial("logit"))
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- glm(limitexp > 1 ~ 
                after + Aword_emperor + Aword_expert + 
                Bword_emperor + Bword_expert, 
              data = dRpanel, family=binomial("logit"))
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt2 <- 
    rbind(c(mbG2$coefficients[2],mbG2c_se[2],
            coefci(mbG2, vcov.=vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id), level=0.95)[2,],
            coefci(mbG2, vcov.=vcovCL(mbG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbG2$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt2 <- as.data.frame(tmpdt2)
  colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt2$frame <- "Biased history"
  tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  ## Combine Data
  mbG1dt <- rbind(tmpdt1,tmpdt2)
  mbG1dt$frame <- factor(mbG1dt$frame, levels = c("Hate speech","Biased history"))
  mbG1dt$endorsement <- factor(mbG1dt$endorsement, levels = unique(mbG1dt$endorsement))
  
  write.csv(mbG1dt, row.names = F, 
            file = paste0("../out/effect_endorsement_logit_base_v2","_",postfixset,".csv"))
  
  ## logit (Extended) ##
  
  tmp1 <- glm(limitexp > 1 ~ 
                after + Aword_expert + Aword_emperor_liberal + 
                Bword_expert + Bword_emperor_liberal + 
                fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
              data = dLpanel, family=binomial("logit"))
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- glm(limitexp ~ 
                after + Aword_emperor + Aword_expert + 
                Bword_emperor + Bword_expert + 
                fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
              data = dLpanel)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt1 <- 
    rbind(c(mG1$coefficients[2],mG1c_se[2],
            coefci(mG1, vcov.=vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id), level=0.95)[2,],
            coefci(mG1, vcov.=vcovCL(mG1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mG1$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt1 <- as.data.frame(tmpdt1)
  colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt1$frame <- "Hate speech"
  tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  tmp1 <- glm(limitexp > 1 ~ 
                after + Aword_expert + Aword_emperor_liberal + 
                Bword_expert + Bword_emperor_liberal + 
                fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
              data = dRpanel, family=binomial("logit"))
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- glm(limitexp ~ 
                after + Aword_emperor + Aword_expert + 
                Bword_emperor + Bword_expert + 
                fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
              data = dRpanel)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt2 <- 
    rbind(c(mG2$coefficients[2],mG2c_se[2],
            coefci(mG2, vcov.=vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id), level=0.95)[2,],
            coefci(mG2, vcov.=vcovCL(mG2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mG2$terms))])$start_id), level=0.90)[2,]),
          c(tmp1$coefficients[2],tmp1c_se[2],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[2,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[2,]),
          c(tmp2$coefficients[2],tmp2c_se[2],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[2,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[2,]))
  
  tmpdt2 <- as.data.frame(tmpdt2)
  colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt2$frame <- "Biased history"
  tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  ## Combine Data
  mG1dt <- rbind(tmpdt1,tmpdt2)
  mG1dt$frame <- factor(mG1dt$frame, levels = c("Hate speech","Biased history"))
  mG1dt$endorsement <- factor(mG1dt$endorsement, levels = unique(mG1dt$endorsement))
  
  write.csv(mG1dt, row.names = F, 
            file = paste0("../out/effect_endorsement_logit_ext_v2","_",postfixset,".csv"))
  
  ## Ordinal logit (Baseline) ##
  
  tmp1 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_expert + Aword_emperor_liberal + 
                 Bword_emperor + Bword_emperor_liberal, 
               data = dLpanel, Hess = TRUE)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_emperor + Aword_expert + 
                 Bword_emperor + Bword_emperor_liberal, 
               data = dLpanel, Hess = TRUE)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt1 <- 
    rbind(c(mbD1$coefficients[1],mbD1c_se[1],
            coefci(mbD1, vcov.=vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id), level=0.95)[1,],
            coefci(mbD1, vcov.=vcovCL(mbD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mbD1$terms))])$start_id), level=0.90)[1,]),
          c(tmp1$coefficients[1],tmp1c_se[1],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
          c(tmp2$coefficients[1],tmp2c_se[1],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))
  
  tmpdt1 <- as.data.frame(tmpdt1)
  colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt1$frame <- "Hate speech"
  tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  tmp1 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_expert + Aword_emperor_liberal + 
                 Bword_emperor + Bword_emperor_liberal, 
               data = dRpanel, Hess = TRUE)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_emperor + Aword_expert + 
                 Bword_emperor + Bword_emperor_liberal, 
               data = dRpanel, Hess = TRUE)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt2 <- 
    rbind(c(mbD2$coefficients[1],mbD2c_se[1],
            coefci(mbD2, vcov.=vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id), level=0.95)[1,],
            coefci(mbD2, vcov.=vcovCL(mbD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mbD2$terms))])$start_id), level=0.90)[1,]),
          c(tmp1$coefficients[1],tmp1c_se[1],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
          c(tmp2$coefficients[1],tmp2c_se[1],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))
  
  tmpdt2 <- as.data.frame(tmpdt2)
  colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt2$frame <- "Biased history"
  tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  ## Combine Data
  mbD1dt <- rbind(tmpdt1,tmpdt2)
  mbD1dt$frame <- factor(mbD1dt$frame, levels = c("Hate speech","Biased history"))
  mbD1dt$endorsement <- factor(mbD1dt$endorsement, levels = unique(mbD1dt$endorsement))
  
  write.csv(mbD1dt, row.names = F, 
            file = paste0("../out/effect_endorsement_ol_base_v2","_",postfixset,".csv"))
  
  ## Ordinal logit (Extended) ##
  
  tmp1 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_expert + Aword_emperor_liberal + 
                 Bword_emperor + Bword_emperor_liberal + 
                 fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
               data = dLpanel, Hess = TRUE)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_emperor + Aword_expert + 
                 Bword_emperor + Bword_emperor_liberal + 
                 fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
               data = dLpanel, Hess = TRUE)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt1 <- 
    rbind(c(mD1$coefficients[1],mD1c_se[1],
            coefci(mD1, vcov.=vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id), level=0.95)[1,],
            coefci(mD1, vcov.=vcovCL(mD1,cluster=na.omit(dLpanel[,c("start_id",all.vars(mD1$terms))])$start_id), level=0.90)[1,]),
          c(tmp1$coefficients[1],tmp1c_se[1],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
          c(tmp2$coefficients[1],tmp2c_se[1],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dLpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))
  
  tmpdt1 <- as.data.frame(tmpdt1)
  colnames(tmpdt1) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt1$frame <- "Hate speech"
  tmpdt1$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  tmp1 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_expert + Aword_emperor_liberal + 
                 Bword_emperor + Bword_emperor_liberal + 
                 fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
               data = dRpanel, Hess = TRUE)
  summary(tmp1)
  coeftest(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id))
  tmp1c_se <- sqrt(diag(vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id)))
  tmp1c_p <- pnorm(-abs(summary(tmp1)$coefficients[,1]/tmp1c_se))*2
  tmp2 <- polr(as.ordered(limitexp) ~ 
                 after + Aword_emperor + Aword_expert + 
                 Bword_emperor + Bword_emperor_liberal + 
                 fem + age + inccat + educat + employed + knall + csup + eqview + hmview, 
               data = dRpanel, Hess = TRUE)
  summary(tmp2)
  coeftest(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id))
  tmp2c_se <- sqrt(diag(vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id)))
  tmp2c_p <- pnorm(-abs(summary(tmp2)$coefficients[,1]/tmp2c_se))*2
  
  tmpdt2 <- 
    rbind(c(mD2$coefficients[1],mD2c_se[1],
            coefci(mD2, vcov.=vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id), level=0.95)[1,],
            coefci(mD2, vcov.=vcovCL(mD2,cluster=na.omit(dRpanel[,c("start_id",all.vars(mD2$terms))])$start_id), level=0.90)[1,]),
          c(tmp1$coefficients[1],tmp1c_se[1],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.95)[1,],
            coefci(tmp1, vcov.=vcovCL(tmp1,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp1$terms))])$start_id), level=0.90)[1,]),
          c(tmp2$coefficients[1],tmp2c_se[1],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.95)[1,],
            coefci(tmp2, vcov.=vcovCL(tmp2,cluster=na.omit(dRpanel[,c("start_id",all.vars(tmp2$terms))])$start_id), level=0.90)[1,]))
  
  tmpdt2 <- as.data.frame(tmpdt2)
  colnames(tmpdt2) <- c("cf","se","lci","uci","lci90","uci90")
  tmpdt2$frame <- "Biased history"
  tmpdt2$endorsement <- c("Experts","Emperor","Liberal\nemperor")
  
  ## Combine Data
  mD1dt <- rbind(tmpdt1,tmpdt2)
  mD1dt$frame <- factor(mD1dt$frame, levels = c("Hate speech","Biased history"))
  mD1dt$endorsement <- factor(mD1dt$endorsement, levels = unique(mD1dt$endorsement))
  
  write.csv(mD1dt, row.names = F, 
            file = paste0("../out/effect_endorsement_ol_ext_v2","_",postfixset,".csv"))
  
  return(NULL)
  
}

## Self-reported
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_self_mod==0)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_self_mod==0)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_self_mod0")
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_self_mod==1)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_self_mod==1)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_self_mod1")

## National Security Ideology
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_iss_1_mod==0)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_iss_1_mod==0)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_iss_1_mod0")
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_iss_1_mod==1)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_iss_1_mod==1)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_iss_1_mod1")

## Equality Ideology
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_iss_2_mod==0)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_iss_2_mod==0)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_iss_2_mod0")
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_iss_2_mod==1)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_iss_2_mod==1)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_iss_2_mod1")

## Party Ideology
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_psup_mod==0)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_psup_mod==0)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_psup_mod0")
dLpaneltmp <- subset(dpanel, frame_left==1 & ide_psup_mod==1)
dRpaneltmp <- subset(dpanel, frame_right==1 & ide_psup_mod==1)
tmpfunc(dLpaneltmp,dRpaneltmp,"ide_psup_mod1")

#'
#' ### Linear Model (Baseline) (Figure 3)
#'

#############
## lm_base ##
#############

dt_selfI <- read.csv("../out/effect_endorsement_lm_base_v2_ide_self_mod0.csv")
dt_selfI$ide <- "Self-reported"
dt_selfI$val <- "Ideologue"
dt_selfM <- read.csv("../out/effect_endorsement_lm_base_v2_ide_self_mod1.csv")
dt_selfM$ide <- "Self-reported"
dt_selfM$val <- "Moderate"
dt_iss_1I <- read.csv("../out/effect_endorsement_lm_base_v2_ide_iss_1_mod0.csv")
dt_iss_1I$ide <- "National\nSecurity"
dt_iss_1I$val <- "Ideologue"
dt_iss_1M <- read.csv("../out/effect_endorsement_lm_base_v2_ide_iss_1_mod1.csv")
dt_iss_1M$ide <- "National\nSecurity"
dt_iss_1M$val <- "Moderate"
dt_iss_2I <- read.csv("../out/effect_endorsement_lm_base_v2_ide_iss_2_mod0.csv")
dt_iss_2I$ide <- "Equality"
dt_iss_2I$val <- "Ideologue"
dt_iss_2M <- read.csv("../out/effect_endorsement_lm_base_v2_ide_iss_2_mod1.csv")
dt_iss_2M$ide <- "Equality"
dt_iss_2M$val <- "Moderate"
dt_psupI <- read.csv("../out/effect_endorsement_lm_base_v2_ide_psup_mod0.csv")
dt_psupI$ide <- "Party"
dt_psupI$val <- "Ideologue"
dt_psupM <- read.csv("../out/effect_endorsement_lm_base_v2_ide_psup_mod1.csv")
dt_psupM$ide <- "Party"
dt_psupM$val <- "Moderate"

dtmod <- rbind(dt_selfI,dt_selfM, 
               dt_iss_1I,dt_iss_1M,
               dt_iss_2I,dt_iss_2M,
               dt_psupI,dt_psupM)
dtmod$frame <- factor(dtmod$frame, levels=unique(dtmod$frame))
dtmod$endorsement <- factor(dtmod$endorsement,levels=unique(dtmod$endorsement))
dtmod$ide <- factor(dtmod$ide, levels=unique(dtmod$ide))
dtmod$val <- factor(dtmod$val, levels=unique(dtmod$val))

p <- ggplot(dtmod, aes(x=ide,y=cf,
                       shape=endorsement,color=endorsement)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.3, position=position_dodge(width=0.5)) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2, position=position_dodge(width=0.5)) + 
  geom_point(size=2.5, position=position_dodge(width=0.5)) + 
  facet_grid(val~frame) + 
  scale_color_manual(name="Endorser", values=rep("black",3)) + 
  scale_shape_discrete(name="Endorser") + 
  labs(x="Type of Ideology Measurement", 
       y="Negative endorsement treatment effect\nOLS coefficient") + 
  theme_bw()

#+ fig.width=8, fig.height=5
p

#+ eval=FALSE
ggsave(paste0("../out/effect_endorsement_lm_base_v2_mod.pdf"), width=8, height=5)
ggsave(paste0("../out/effect_endorsement_lm_base_v2_mod.png"), width=8, height=5)

#'
#' ### Linear Model (Extended) (Figure H.1)
#'

############
## lm_ext ##
############

dt_selfI <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_self_mod0.csv")
dt_selfI$ide <- "Self-reported"
dt_selfI$val <- "Ideologue"
dt_selfM <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_self_mod1.csv")
dt_selfM$ide <- "Self-reported"
dt_selfM$val <- "Moderate"
dt_iss_1I <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_iss_1_mod0.csv")
dt_iss_1I$ide <- "National\nSecurity"
dt_iss_1I$val <- "Ideologue"
dt_iss_1M <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_iss_1_mod1.csv")
dt_iss_1M$ide <- "National\nSecurity"
dt_iss_1M$val <- "Moderate"
dt_iss_2I <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_iss_2_mod0.csv")
dt_iss_2I$ide <- "Equality"
dt_iss_2I$val <- "Ideologue"
dt_iss_2M <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_iss_2_mod1.csv")
dt_iss_2M$ide <- "Equality"
dt_iss_2M$val <- "Moderate"
dt_psupI <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_psup_mod0.csv")
dt_psupI$ide <- "Party"
dt_psupI$val <- "Ideologue"
dt_psupM <- read.csv("../out/effect_endorsement_lm_ext_v2_ide_psup_mod1.csv")
dt_psupM$ide <- "Party"
dt_psupM$val <- "Moderate"

dtmod <- rbind(dt_selfI,dt_selfM, 
               dt_iss_1I,dt_iss_1M,
               dt_iss_2I,dt_iss_2M,
               dt_psupI,dt_psupM)
dtmod$frame <- factor(dtmod$frame, levels=unique(dtmod$frame))
dtmod$endorsement <- factor(dtmod$endorsement,levels=unique(dtmod$endorsement))
dtmod$ide <- factor(dtmod$ide, levels=unique(dtmod$ide))
dtmod$val <- factor(dtmod$val, levels=unique(dtmod$val))

p <- ggplot(dtmod, aes(x=ide,y=cf,
                       shape=endorsement,color=endorsement)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.3, position=position_dodge(width=0.5)) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2, position=position_dodge(width=0.5)) + 
  geom_point(size=2.5, position=position_dodge(width=0.5)) + 
  facet_grid(val~frame) + 
  scale_color_manual(name="Endorser", values=rep("black",3)) + 
  scale_shape_discrete(name="Endorser") + 
  labs(x="Type of Ideology Measurement", 
       y="Negative endorsement treatment effect\nOLS coefficient") + 
  theme_bw()

#+ fig.width=8, fig.height=5
p

#+ eval=FALSE
ggsave(paste0("../out/effect_endorsement_lm_ext_v2_mod.pdf"), width=8, height=5)
ggsave(paste0("../out/effect_endorsement_lm_ext_v2_mod.png"), width=8, height=5)

#'
#' ### Logit (Baseline) (Figure H.2)
#'

################
## logit_base ##
################

dt_selfI <- read.csv("../out/effect_endorsement_logit_base_v2_ide_self_mod0.csv")
dt_selfI$ide <- "Self-reported"
dt_selfI$val <- "Ideologue"
dt_selfM <- read.csv("../out/effect_endorsement_logit_base_v2_ide_self_mod1.csv")
dt_selfM$ide <- "Self-reported"
dt_selfM$val <- "Moderate"
dt_iss_1I <- read.csv("../out/effect_endorsement_logit_base_v2_ide_iss_1_mod0.csv")
dt_iss_1I$ide <- "National\nSecurity"
dt_iss_1I$val <- "Ideologue"
dt_iss_1M <- read.csv("../out/effect_endorsement_logit_base_v2_ide_iss_1_mod1.csv")
dt_iss_1M$ide <- "National\nSecurity"
dt_iss_1M$val <- "Moderate"
dt_iss_2I <- read.csv("../out/effect_endorsement_logit_base_v2_ide_iss_2_mod0.csv")
dt_iss_2I$ide <- "Equality"
dt_iss_2I$val <- "Ideologue"
dt_iss_2M <- read.csv("../out/effect_endorsement_logit_base_v2_ide_iss_2_mod1.csv")
dt_iss_2M$ide <- "Equality"
dt_iss_2M$val <- "Moderate"
dt_psupI <- read.csv("../out/effect_endorsement_logit_base_v2_ide_psup_mod0.csv")
dt_psupI$ide <- "Party"
dt_psupI$val <- "Ideologue"
dt_psupM <- read.csv("../out/effect_endorsement_logit_base_v2_ide_psup_mod1.csv")
dt_psupM$ide <- "Party"
dt_psupM$val <- "Moderate"

dtmod <- rbind(dt_selfI,dt_selfM, 
               dt_iss_1I,dt_iss_1M,
               dt_iss_2I,dt_iss_2M,
               dt_psupI,dt_psupM)
dtmod$frame <- factor(dtmod$frame, levels=unique(dtmod$frame))
dtmod$endorsement <- factor(dtmod$endorsement,levels=unique(dtmod$endorsement))
dtmod$ide <- factor(dtmod$ide, levels=unique(dtmod$ide))
dtmod$val <- factor(dtmod$val, levels=unique(dtmod$val))

p <- ggplot(dtmod, aes(x=ide,y=cf,
                       shape=endorsement,color=endorsement)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.3, position=position_dodge(width=0.5)) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2, position=position_dodge(width=0.5)) + 
  geom_point(size=2.5, position=position_dodge(width=0.5)) + 
  facet_grid(val~frame) + 
  scale_color_manual(name="Endorser", values=rep("black",3)) + 
  scale_shape_discrete(name="Endorser") + 
  labs(x="Type of Ideology Measurement", 
       y="Negative endorsement treatment effect\nlogit coefficient") + 
  theme_bw()

#+ fig.width=8, fig.height=5
p

#+ eval=FALSE
ggsave(paste0("../out/effect_endorsement_logit_base_v2_mod.pdf"), width=8, height=5)
ggsave(paste0("../out/effect_endorsement_logit_base_v2_mod.png"), width=8, height=5)

#'
#' ### Logit (Extended) (Figure H.3)
#'

###############
## logit_ext ##
###############

dt_selfI <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_self_mod0.csv")
dt_selfI$ide <- "Self-reported"
dt_selfI$val <- "Ideologue"
dt_selfM <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_self_mod1.csv")
dt_selfM$ide <- "Self-reported"
dt_selfM$val <- "Moderate"
dt_iss_1I <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_iss_1_mod0.csv")
dt_iss_1I$ide <- "National\nSecurity"
dt_iss_1I$val <- "Ideologue"
dt_iss_1M <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_iss_1_mod1.csv")
dt_iss_1M$ide <- "National\nSecurity"
dt_iss_1M$val <- "Moderate"
dt_iss_2I <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_iss_2_mod0.csv")
dt_iss_2I$ide <- "Equality"
dt_iss_2I$val <- "Ideologue"
dt_iss_2M <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_iss_2_mod1.csv")
dt_iss_2M$ide <- "Equality"
dt_iss_2M$val <- "Moderate"
dt_psupI <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_psup_mod0.csv")
dt_psupI$ide <- "Party"
dt_psupI$val <- "Ideologue"
dt_psupM <- read.csv("../out/effect_endorsement_logit_ext_v2_ide_psup_mod1.csv")
dt_psupM$ide <- "Party"
dt_psupM$val <- "Moderate"

dtmod <- rbind(dt_selfI,dt_selfM, 
               dt_iss_1I,dt_iss_1M,
               dt_iss_2I,dt_iss_2M,
               dt_psupI,dt_psupM)
dtmod$frame <- factor(dtmod$frame, levels=unique(dtmod$frame))
dtmod$endorsement <- factor(dtmod$endorsement,levels=unique(dtmod$endorsement))
dtmod$ide <- factor(dtmod$ide, levels=unique(dtmod$ide))
dtmod$val <- factor(dtmod$val, levels=unique(dtmod$val))

p <- ggplot(dtmod, aes(x=ide,y=cf,
                       shape=endorsement,color=endorsement)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.3, position=position_dodge(width=0.5)) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2, position=position_dodge(width=0.5)) + 
  geom_point(size=2.5, position=position_dodge(width=0.5)) + 
  facet_grid(val~frame) + 
  scale_color_manual(name="Endorser", values=rep("black",3)) + 
  scale_shape_discrete(name="Endorser") + 
  labs(x="Type of Ideology Measurement", 
       y="Negative endorsement treatment effect\nlogit coefficient") + 
  theme_bw()

#+ fig.width=8, fig.height=5
p

#+ eval=FALSE
ggsave(paste0("../out/effect_endorsement_logit_ext_v2_mod.pdf"), width=8, height=5)
ggsave(paste0("../out/effect_endorsement_logit_ext_v2_mod.png"), width=8, height=5)

#'
#' ### Ordinal Logit (Baseline) (Figure H.4)
#'

#############
## ol_base ##
#############

dt_selfI <- read.csv("../out/effect_endorsement_ol_base_v2_ide_self_mod0.csv")
dt_selfI$ide <- "Self-reported"
dt_selfI$val <- "Ideologue"
dt_selfM <- read.csv("../out/effect_endorsement_ol_base_v2_ide_self_mod1.csv")
dt_selfM$ide <- "Self-reported"
dt_selfM$val <- "Moderate"
dt_iss_1I <- read.csv("../out/effect_endorsement_ol_base_v2_ide_iss_1_mod0.csv")
dt_iss_1I$ide <- "National\nSecurity"
dt_iss_1I$val <- "Ideologue"
dt_iss_1M <- read.csv("../out/effect_endorsement_ol_base_v2_ide_iss_1_mod1.csv")
dt_iss_1M$ide <- "National\nSecurity"
dt_iss_1M$val <- "Moderate"
dt_iss_2I <- read.csv("../out/effect_endorsement_ol_base_v2_ide_iss_2_mod0.csv")
dt_iss_2I$ide <- "Equality"
dt_iss_2I$val <- "Ideologue"
dt_iss_2M <- read.csv("../out/effect_endorsement_ol_base_v2_ide_iss_2_mod1.csv")
dt_iss_2M$ide <- "Equality"
dt_iss_2M$val <- "Moderate"
dt_psupI <- read.csv("../out/effect_endorsement_ol_base_v2_ide_psup_mod0.csv")
dt_psupI$ide <- "Party"
dt_psupI$val <- "Ideologue"
dt_psupM <- read.csv("../out/effect_endorsement_ol_base_v2_ide_psup_mod1.csv")
dt_psupM$ide <- "Party"
dt_psupM$val <- "Moderate"

dtmod <- rbind(dt_selfI,dt_selfM, 
               dt_iss_1I,dt_iss_1M,
               dt_iss_2I,dt_iss_2M,
               dt_psupI,dt_psupM)
dtmod$frame <- factor(dtmod$frame, levels=unique(dtmod$frame))
dtmod$endorsement <- factor(dtmod$endorsement,levels=unique(dtmod$endorsement))
dtmod$ide <- factor(dtmod$ide, levels=unique(dtmod$ide))
dtmod$val <- factor(dtmod$val, levels=unique(dtmod$val))

p <- ggplot(dtmod, aes(x=ide,y=cf,
                       shape=endorsement,color=endorsement)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.3, position=position_dodge(width=0.5)) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2, position=position_dodge(width=0.5)) + 
  geom_point(size=2.5, position=position_dodge(width=0.5)) + 
  facet_grid(val~frame) + 
  scale_color_manual(name="Endorser", values=rep("black",3)) + 
  scale_shape_discrete(name="Endorser") + 
  labs(x="Type of Ideology Measurement", 
       y="Negative endorsement treatment effect\nordinal logit coefficient") + 
  theme_bw()

#+ fig.width=8, fig.height=5
p

#+ eval=FALSE
ggsave(paste0("../out/effect_endorsement_ol_base_v2_mod.pdf"), width=8, height=5)
ggsave(paste0("../out/effect_endorsement_ol_base_v2_mod.png"), width=8, height=5)

#'
#' ### Ordinal Logit (Extended) (Figure H.5)
#'

############
## ol_ext ##
############

dt_selfI <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_self_mod0.csv")
dt_selfI$ide <- "Self-reported"
dt_selfI$val <- "Ideologue"
dt_selfM <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_self_mod1.csv")
dt_selfM$ide <- "Self-reported"
dt_selfM$val <- "Moderate"
dt_iss_1I <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_iss_1_mod0.csv")
dt_iss_1I$ide <- "National\nSecurity"
dt_iss_1I$val <- "Ideologue"
dt_iss_1M <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_iss_1_mod1.csv")
dt_iss_1M$ide <- "National\nSecurity"
dt_iss_1M$val <- "Moderate"
dt_iss_2I <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_iss_2_mod0.csv")
dt_iss_2I$ide <- "Equality"
dt_iss_2I$val <- "Ideologue"
dt_iss_2M <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_iss_2_mod1.csv")
dt_iss_2M$ide <- "Equality"
dt_iss_2M$val <- "Moderate"
dt_psupI <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_psup_mod0.csv")
dt_psupI$ide <- "Party"
dt_psupI$val <- "Ideologue"
dt_psupM <- read.csv("../out/effect_endorsement_ol_ext_v2_ide_psup_mod1.csv")
dt_psupM$ide <- "Party"
dt_psupM$val <- "Moderate"

dtmod <- rbind(dt_selfI,dt_selfM, 
               dt_iss_1I,dt_iss_1M,
               dt_iss_2I,dt_iss_2M,
               dt_psupI,dt_psupM)
dtmod$frame <- factor(dtmod$frame, levels=unique(dtmod$frame))
dtmod$endorsement <- factor(dtmod$endorsement,levels=unique(dtmod$endorsement))
dtmod$ide <- factor(dtmod$ide, levels=unique(dtmod$ide))
dtmod$val <- factor(dtmod$val, levels=unique(dtmod$val))

p <- ggplot(dtmod, aes(x=ide,y=cf,
                       shape=endorsement,color=endorsement)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lci,ymax=uci), width=0.3, position=position_dodge(width=0.5)) + 
  geom_errorbar(aes(ymin=lci90,ymax=uci90), width=0, size=1.2, position=position_dodge(width=0.5)) + 
  geom_point(size=2.5, position=position_dodge(width=0.5)) + 
  facet_grid(val~frame) + 
  scale_color_manual(name="Endorser", values=rep("black",3)) + 
  scale_shape_discrete(name="Endorser") + 
  labs(x="Type of Ideology Measurement", 
       y="Negative endorsement treatment effect\nordinal logit coefficient") + 
  theme_bw()

#+ fig.width=8, fig.height=5
p

#+ eval=FALSE
ggsave(paste0("../out/effect_endorsement_ol_ext_v2_mod.pdf"), width=8, height=5)
ggsave(paste0("../out/effect_endorsement_ol_ext_v2_mod.png"), width=8, height=5)

#+ eval=FALSE, echo=FALSE
# Exporting PDF File
# In R Studio
# rmarkdown::render('main_analysis_emperor_final.R', 'pdf_document', encoding = 'UTF-8')
# In Terminal, run:
# Rscript -e "rmarkdown::render('main_analysis_emperor_final.R', 'pdf_document', encoding = 'UTF-8')"


